diff --git a/code/IO.f90 b/code/IO.f90 index 460831e98..0e768bf02 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -168,14 +168,14 @@ end function if (FEsolver == 'Abaqus') then open(unit+1,status='old',err=100,& file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//InputFileExtension) + trim(getModelName())//InputFileExtension) open(unit,err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//InputFileExtension//'_assembly') + trim(getModelName())//InputFileExtension//'_assembly') IO_open_inputFile = IO_abaqus_assembleInputFile(unit,unit+1) ! strip comments and concatenate any "include"s close(unit+1) else open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//InputFileExtension) + trim(getModelName())//InputFileExtension) IO_open_inputFile = .true. endif @@ -1142,10 +1142,12 @@ endfunction case (45) msg = 'error opening spectral loadcase' case (46) - msg = 'missing parameter in spectral loadcase' - case (47) msg = 'mask consistency violated in spectral loadcase' + case (47) + msg = 'negative time increment in spectral loadcase' case (48) + msg = 'Non-positive increment in spectral loadcase' + case (49) msg = 'Non-positive relative parameter for spectral method' case (50) msg = 'Error writing constitutive output description' @@ -1155,6 +1157,10 @@ endfunction msg = 'Cannot open input file' case (102) msg = 'argument count error (mesh and loadcase) for mpie_spectral' + case (103) + msg = 'Resolution contains odd number' + case (104) + msg = 'FFTW init error' case (105) msg = 'Error reading from ODF file' case (110) @@ -1308,6 +1314,9 @@ endfunction case (700) msg = 'Singular matrix in stress iteration' + case (800) + msg = 'Matrix inversion error' + ! Error messages related to parsing of Abaqus input file case (900) msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)' diff --git a/code/libfftw3.a b/code/libfftw3.a new file mode 100644 index 000000000..80a58dadf Binary files /dev/null and b/code/libfftw3.a differ diff --git a/code/libfftw3_threads.a b/code/libfftw3_threads.a new file mode 100644 index 000000000..29ac934cc Binary files /dev/null and b/code/libfftw3_threads.a differ diff --git a/code/makefile b/code/makefile index d12a973b9..f8ce20266 100644 --- a/code/makefile +++ b/code/makefile @@ -6,11 +6,13 @@ # Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads" and "make", "make install" is not needed # as long as the two library files are copied to the source code directory. -cpspectral.out: mpie_spectral.o CPFEM.a - ifort -openmp -o cpspectral.out mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread +cpspectral.exe: mpie_spectral.o CPFEM.a + ifort -openmp -o cpspectral.exe mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread + rm *.mod + mpie_spectral.o: mpie_spectral.f90 CPFEM.o ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90 - + CPFEM.a: CPFEM.o ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o diff --git a/code/mpie_cpfem_abaqus_exp.f b/code/mpie_cpfem_abaqus_exp.f index b0155ee05..2e3ec4d21 100644 --- a/code/mpie_cpfem_abaqus_exp.f +++ b/code/mpie_cpfem_abaqus_exp.f @@ -24,7 +24,9 @@ character(len=4), parameter :: InputFileExtension = '.inp' CONTAINS +!-------------------- subroutine mpie_interface_init() +!-------------------- write(6,*) write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>' write(6,*) '$Id$' @@ -32,7 +34,9 @@ subroutine mpie_interface_init() return end subroutine +!-------------------- function getSolverWorkingDirectoryName() +!-------------------- use prec implicit none character(1024) getSolverWorkingDirectoryName @@ -40,10 +44,21 @@ function getSolverWorkingDirectoryName() getSolverWorkingDirectoryName='' CALL VGETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR ) -! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName end function + +!-------------------- +function getModelName() +!-------------------- + character(1024) getModelName + + getModelName = getSolverJobName() +end function + + +!-------------------- function getSolverJobName() +!-------------------- use prec implicit none @@ -52,7 +67,6 @@ function getSolverJobName() getSolverJobName='' CALL VGETJOBNAME(getSolverJobName , LENJOBNAME ) -! write(6,*) 'getSolverJobName', getSolverJobName end function END MODULE diff --git a/code/mpie_cpfem_abaqus_std.f b/code/mpie_cpfem_abaqus_std.f index 7361c284a..1fd8f633f 100644 --- a/code/mpie_cpfem_abaqus_std.f +++ b/code/mpie_cpfem_abaqus_std.f @@ -24,7 +24,9 @@ character(len=4), parameter :: InputFileExtension = '.inp' CONTAINS +!-------------------- subroutine mpie_interface_init() +!-------------------- write(6,*) write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>' write(6,*) '$Id$' @@ -32,7 +34,9 @@ subroutine mpie_interface_init() return end subroutine +!-------------------- function getSolverWorkingDirectoryName() +!-------------------- use prec implicit none character(1024) getSolverWorkingDirectoryName @@ -43,7 +47,19 @@ function getSolverWorkingDirectoryName() ! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName end function + +!-------------------- +function getModelName() +!-------------------- + character(1024) getModelName + + getModelName = getSolverJobName() +end function + + +!-------------------- function getSolverJobName() +!-------------------- use prec implicit none diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index f76044c73..af56dfd9d 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -71,6 +71,17 @@ function getSolverWorkingDirectoryName() ! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName end function + +function getModelName() + use prec + implicit none + + character(1024) getModelName + + getModelName = getSolverJobName() +end function + + function getSolverJobName() use prec implicit none diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 513eeb9ff..6e3017bef 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -30,7 +30,7 @@ program mpie_spectral use math use CPFEM, only: CPFEM_general, CPFEM_initAll use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol,& - itmax, fast_execution, mpieNumThreadsInt + itmax, memory_efficient, mpieNumThreadsInt use homogenization, only: materialpoint_sizeResults, materialpoint_results !$ use OMP_LIB ! the openMP function library @@ -78,7 +78,7 @@ program mpie_spectral complex(pReal), dimension(3,3) :: temp33_Complex real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat - real(pReal), dimension(3) :: xi, xi_middle + real(pReal), dimension(3) :: xi, xi_central integer(pInt), dimension(3) :: k_s integer*8, dimension(2) :: plan_fft @@ -109,8 +109,9 @@ program mpie_spectral ! Reading the loadcase file and assign variables path = getLoadcaseName() - print*,'Loadcase: ',trim(path) - print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) + print '(a,/,a)', 'Loadcase: ',trim(path) + print '(a,/,a)', 'Workingdir: ',trim(getSolverWorkingDirectoryName()) + print '(a,/,a)', 'SolverJobName: ',trim(getSolverJobName()) if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg = path) @@ -131,8 +132,6 @@ program mpie_spectral N_n = N_n+1 end select enddo ! count all identifiers to allocate memory and do sanity check - if ((N_l /= N_s).or.(N_s /= N_t).or.(N_t /= N_n)) & ! sanity check - call IO_error(46,ext_msg = path) ! error message for incomplete input file enddo 101 N_Loadcases = N_l @@ -150,7 +149,7 @@ program mpie_spectral read(unit,'(a1024)',END = 200) line if (IO_isBlank(line)) cycle ! skip empty lines i = i + 1 - posInput = IO_stringPos(line,maxNchunksInput) ! ToDo: Add error message for case that information is not complete + posInput = IO_stringPos(line,maxNchunksInput) do j = 1,maxNchunksInput,2 select case (IO_lc(IO_stringValue(line,posInput,j))) case('l','velocitygrad') @@ -179,20 +178,22 @@ program mpie_spectral 200 close(unit) do i = 1, N_Loadcases - if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(47,i) ! bc_mask consistency - print '(a,/,3(3(f12.6,x)/))','L' ,math_transpose3x3(bc_velocityGrad(:,:,i)) - print '(a,/,3(3(f12.6,x)/))','bc_stress',math_transpose3x3(bc_stress(:,:,i)) - print '(a,/,3(3(l,x)/))', 'bc_mask for velocitygrad',transpose(bc_mask(:,:,1,i)) - print '(a,/,3(3(l,x)/))', 'bc_mask for stress' ,transpose(bc_mask(:,:,2,i)) - print *,'time',bc_timeIncrement(i) - print *,'incs',bc_steps(i) + if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(46,i) ! bc_mask consistency + if (bc_timeIncrement(i) < 0.0_pReal) call IO_error(47,i) ! negative time increment forbidden + if (bc_steps(i) < 1_pInt) call IO_error(48,i) ! non-positive increments requested + print '(a,/,3(3(f12.6,x)/))','L:' ,math_transpose3x3(bc_velocityGrad(:,:,i)) + print '(a,/,3(3(f12.6,x)/))','bc_stress:',math_transpose3x3(bc_stress(:,:,i)) + print '(a,/,3(3(l,x)/))', 'bc_mask for velocitygrad:',transpose(bc_mask(:,:,1,i)) + print '(a,/,3(3(l,x)/))', 'bc_mask for stress:' ,transpose(bc_mask(:,:,2,i)) + print '(a,f12.6)','time: ',bc_timeIncrement(i) + print '(a,i5)','incs: ',bc_steps(i) print *, '' enddo !read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90 - path = getSolverJobName() - print*,'JobName: ',trim(path) - if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = path) + path = getModelName() + print '(a,a)', 'GeomName: ',trim(path) + if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = trim(path)//InputFileExtension) rewind(unit) do @@ -233,11 +234,11 @@ program mpie_spectral enddo 100 close(unit) - if(mod(resolution(1),2)/=0 .or. mod(resolution(2),2)/=0 .or. mod(resolution(3),2)/=0) call IO_error(102) !!ToDo: add correct error to IO + if(mod(resolution(1),2)/=0 .or. mod(resolution(2),2)/=0 .or. mod(resolution(3),2)/=0) call IO_error(103) - print '(a,/,i4,i4,i4)','resolution a b c', resolution - print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension - print *,'homogenization',homog + print '(a,/,i4,i4,i4)','resolution a b c:', resolution + print '(a,/,f6.1,f6.1,f6.1)','dimension x y z:', geomdimension + print '(a,i4)','homogenization: ',homog allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal @@ -261,10 +262,12 @@ program mpie_spectral c066 = c066 * wgt c0 = math_mandel66to3333(c066) call math_invert(6, c066, s066,i, errmatinv) + if(errmatinv) call IO_error(800) !Matrix inversion error s0 = math_mandel66to3333(s066) -!calculation of calculate gamma_hat field in the case of fast execution (needs a lot of memory) - if(fast_execution) then + if(memory_efficient) then ! allocate just single fourth order tensor + allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal + else ! precalculation of gamma_hat field allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal do k = 1, resolution(3) k_s(3) = k-1 @@ -280,32 +283,31 @@ program mpie_spectral xi(1) = real(k_s(1), pReal)/geomdimension(1) if (any(xi /= 0.0_pReal)) then do l = 1,3; do m = 1,3 - xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2) + xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2) enddo; enddo + temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad)) else - xinormdyad = 0.0_pReal + xinormdyad = 0.0_pReal + temp33_Real = 0.0_pReal endif - temp33_Real = math_mul3333xx33(c0, xinormdyad) - temp33_Real = math_inv3x3(temp33_Real) do l=1,3; do m=1,3; do n=1,3; do p=1,3 - gamma_hat(i,j,k, l,m,n,p) = - (0.5*temp33_Real(l,n)+0.5*temp33_Real(n,l)) *& - (0.5*xinormdyad(m,p)+0.5*xinormdyad(p,m)) + gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *& + (xinormdyad(m,p)+xinormdyad(p,m)) enddo; enddo; enddo; enddo enddo; enddo; enddo - else ! or allocate just one fourth order tensor - allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal endif -! calculate xi for the calculation of divergence in Fourier space (middle frequency) - xi_middle(3) = 0.0_pReal !2D case - if(resolution(3) > 1) xi_middle(3) = real(resolution(3)/2, pReal)/geomdimension(3) !3D case - xi_middle(2) = real(resolution(2)/2, pReal)/geomdimension(2) - xi_middle(1) = real(resolution(1)/2, pReal)/geomdimension(1) +! calculate xi for the calculation of divergence in Fourier space (central frequency) + xi_central(3) = 0.0_pReal !2D case + if(resolution(3) > 1) xi_central(3) = real(resolution(3)/2, pReal)/geomdimension(3) !3D case + xi_central(2) = real(resolution(2)/2, pReal)/geomdimension(2) + xi_central(1) = real(resolution(1)/2, pReal)/geomdimension(1) allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal ! Initialization of fftw (see manual on fftw.org for more details) - call dfftw_init_threads(ierr) !toDo: add error code + call dfftw_init_threads(ierr) + if(ierr == 0) call IO_error(104,ierr) call dfftw_plan_with_nthreads(mpieNumThreadsInt) call dfftw_plan_many_dft_r2c(plan_fft(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),& @@ -316,7 +318,6 @@ program mpie_spectral ! write header of output file open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& - //'_'//trim(getLoadcase())& //'.spectralOut',form='UNFORMATTED') write(538), 'load',trim(getLoadcaseName()) write(538), 'workingdir',trim(getSolverWorkingDirectoryName()) @@ -436,10 +437,10 @@ program mpie_spectral enddo; enddo err_div = 2 * err_div_tol err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) - print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av) - print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ' ,math_transpose3x3(cstress_av)/1.e6 - print '(2(a,E8.2))', ' error stress ',err_stress, ' Tol. = ', err_stress_tol - print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8 + print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient:',math_transpose3x3(defgrad_av) + print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress / MPa: ' ,math_transpose3x3(cstress_av)/1.e6 + print '(2(a,E8.2))', ' error stress: ',err_stress, ' Tol. = ', err_stress_tol + print '(2(a,E8.2))', ' error deformation gradient: ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8 if(err_stress < err_stress_tol*0.8) then calcmode = 1_pInt endif @@ -457,10 +458,11 @@ program mpie_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt - call CPFEM_general(2,& + call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) + CPFEM_mode = 2_pInt workfft(i,j,k,:,:) = pstress cstress_av = cstress_av + math_mandel6to33(cstress) enddo; enddo; enddo @@ -476,20 +478,10 @@ program mpie_spectral sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:) + (workfft(2,1,1,m,:))*img))) enddo err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)+1,resolution(2)/2+1,resolution(3)/2+1,:,:)+& ! L infinity Norm of div(stress) - workfft(resolution(1)+2,resolution(2)/2+1,resolution(3)/2+1,:,:)*img,xi_middle)))) + workfft(resolution(1)+2,resolution(2)/2+1,resolution(3)/2+1,:,:)*img,xi_central)))) err_div = err_div/sigma0 !weighting of error - if(fast_execution) then ! fast execution with stored gamma_hat - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 - do m = 1,3; do n = 1,3 - temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)& - + workfft(i*2 ,j,k,:,:)*img)) - enddo; enddo - workfft(i*2-1,j,k,:,:) = real (temp33_Complex) - workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) - enddo; enddo; enddo - - else ! memory saving version, in-time calculation of gamma_hat + if(memory_efficient) then ! memory saving version, in-time calculation of gamma_hat do k = 1, resolution(3) k_s(3) = k-1 if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) @@ -501,19 +493,19 @@ program mpie_spectral xi(3) = 0.0_pReal !for the 2D case if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) !3D case xi(2) = real(k_s(2), pReal)/geomdimension(2) - xi(1) = real(k_s(1), pReal)/geomdimension(1) + xi(1) = real(k_s(1), pReal)/geomdimension(1) if (any(xi(:) /= 0.0_pReal)) then do l = 1,3; do m = 1,3 - xinormdyad(l,m) = xi(l)*xi( m)/sum(xi**2) + xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2) enddo; enddo + temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad)) else xinormdyad = 0.0_pReal + temp33_Real = 0.0_pReal endif - temp33_Real = math_mul3333xx33(c0, xinormdyad) - temp33_Real = math_inv3x3(temp33_Real) do l=1,3; do m=1,3; do n=1,3; do p=1,3 - gamma_hat(1,1,1, l,m,n,p) = - (0.5*temp33_Real(l,n)+0.5*temp33_Real(n,l))*& - (0.5*xinormdyad(m,p) +0.5*xinormdyad(p,m)) + gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*& + (xinormdyad(m,p) +xinormdyad(p,m)) enddo; enddo; enddo; enddo do m = 1,3; do n = 1,3 temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,:,:) *(workfft(i*2-1,j,k,:,:)& @@ -522,6 +514,15 @@ program mpie_spectral workfft(i*2-1,j,k,:,:) = real (temp33_Complex) workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) enddo; enddo; enddo + else !use precalculated gamma-operator + do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 + do m = 1,3; do n = 1,3 + temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)& + + workfft(i*2 ,j,k,:,:)*img)) + enddo; enddo + workfft(i*2-1,j,k,:,:) = real (temp33_Complex) + workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) + enddo; enddo; enddo endif workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency (real part) @@ -538,9 +539,9 @@ program mpie_spectral err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel !accecpt relativ error specified err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) - print '(2(a,E8.2))', ' error divergence ',err_div, ' Tol. = ', err_div_tol - print '(2(a,E8.2))', ' error stress ',err_stress, ' Tol. = ', err_stress_tol - print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol + print '(2(a,E8.2))', ' error divergence: ',err_div, ' Tol. = ', err_div_tol + print '(2(a,E8.2))', ' error stress: ',err_stress, ' Tol. = ', err_stress_tol + print '(2(a,E8.2))', ' error deformation gradient: ',err_defgrad,' Tol. = ', err_defgrad_tol if((err_stress > err_stress_tol .or. err_defgrad > err_defgrad_tol) .and. err_div < err_div_tol) then ! change to calculation of BCs, reset damper etc. calcmode = 0_pInt @@ -552,10 +553,10 @@ program mpie_spectral write(538) materialpoint_results(:,1,:) !write to output file - print '(a,x,f12.7)' , ' Determinant of Deformation Aim:', math_det3x3(defgradAim) - print '(a,/,3(3(f12.7,x)/))', ' Deformation Aim: ',math_transpose3x3(defgradAim) - print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av) - print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',math_transpose3x3(cstress_av)/1.e6 + print '(a,x,f12.7)' , ' Determinant of Deformation Aim: ', math_det3x3(defgradAim) + print '(a,/,3(3(f12.7,x)/))', ' Deformation Aim: ',math_transpose3x3(defgradAim) + print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient:',math_transpose3x3(defgrad_av) + print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress / MPa: ',math_transpose3x3(cstress_av)/1.e6 print '(A)', '************************************************************' enddo ! end looping over steps in current loadcase enddo ! end looping over loadcases diff --git a/code/mpie_spectral_interface.f90 b/code/mpie_spectral_interface.f90 index 4a4c57801..c7703ab09 100644 --- a/code/mpie_spectral_interface.f90 +++ b/code/mpie_spectral_interface.f90 @@ -59,28 +59,43 @@ function getSolverJobName() implicit none - character(1024) getSolverJobName, outName, cwd + character(1024) getSolverJobName + getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase()) + +endfunction + +!******************************************************************** +! basename of geometry file from command line arguments +! +!******************************************************************** +function getModelName() + + use prec, only: pInt + + implicit none + + character(1024) getModelName, outName, cwd character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ integer(pInt) posExt,posSep - getSolverJobName = '' + getModelName = '' call getarg(1,outName) posExt = scan(outName,'.',back=.true.) posSep = scan(outName,pathSep,back=.true.) if (posExt <= posSep) posExt = len_trim(outName)+1 ! no extension present - getSolverJobName = outName(1:posExt-1) ! path to geometry file (excl. extension) + getModelName = outName(1:posExt-1) ! path to geometry file (excl. extension) - if (scan(getSolverJobName,pathSep) /= 1) then ! relative path given as command line argument + if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument call getcwd(cwd) - getSolverJobName = rectifyPath(trim(cwd)//'/'//getSolverJobName) + getModelName = rectifyPath(trim(cwd)//'/'//getModelName) else - getSolverJobName = rectifyPath(getSolverJobName) + getModelName = rectifyPath(getModelName) endif - getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& - getSolverJobName) + getModelName = makeRelativePath(getSolverWorkingDirectoryName(),& + getModelName) return endfunction diff --git a/code/numerics.f90 b/code/numerics.f90 index 0c43dacab..8172332ba 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -50,7 +50,7 @@ real(pReal) relevantStrain, & ! strain err_stress_tol, & ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value) err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol err_defgrad_tol ! tolerance for error of defgrad compared to prescribed defgrad -logical fast_execution ! for fast execution (pre calculation of gamma_hat) +logical memory_efficient ! for fast execution (pre calculation of gamma_hat) integer(pInt) itmax , & ! maximum number of iterations @@ -94,7 +94,7 @@ subroutine numerics_init() character(len=1024) line ! OpenMP variable -!$ character(len=2) mpieNumThreadsString !enironment variable MPIE_NUMTHREADS +!$ character(len=4) mpieNumThreadsString !environment variable MPIE_NUMTHREADS write(6,*) write(6,*) '<<<+- numerics init -+>>>' @@ -142,11 +142,11 @@ subroutine numerics_init() volDiscrPow_RGC = 5.0 !* spectral parameters: - err_div_tol = 1.0e-4 - err_defgrad_tol = 1.0e-3 - err_stress_tolrel = 0.01 - itmax = 20_pInt - fast_execution = .false. + err_div_tol = 1.0e-4 ! proposed by Suquet, less strict criteria are usefull, e.g. 5e-3 + err_defgrad_tol = 1.0e-3 ! relative tolerance for fullfillment of average deformation gradient (is usually passively fullfilled) + err_stress_tolrel = 0.01 ! relative tolerance for fullfillment of stress BC + itmax = 20_pInt ! Maximum iteration number + memory_efficient = .true. ! Precalculate Gamma-operator (81 double per point) !* Random seeding parameters: added <<>> fixedSeed = 0_pInt @@ -154,7 +154,7 @@ subroutine numerics_init() !* determin number of threads from environment variable MPIE_NUM_THREADS !$ call GetEnv('MPIE_NUM_THREADS',mpieNumThreadsString) ! get environment variable MPIE_NUM_THREADS... -!$ read(mpieNumThreadsString,'(i2)') mpieNumThreadsInt ! ...convert it to integer... +!$ read(mpieNumThreadsString,'(i4)') mpieNumThreadsInt ! ...convert it to integer... !$ if (mpieNumThreadsInt < 1) mpieNumThreadsInt = 1 ! ...ensure that its at least one... !$ call omp_set_num_threads(mpieNumThreadsInt) ! ...and use it as number of threads for parallel execution @@ -256,8 +256,8 @@ subroutine numerics_init() err_stress_tolrel = IO_floatValue(line,positions,2) case ('itmax') itmax = IO_intValue(line,positions,2) - case ('fast_execution') - fast_execution = IO_intValue(line,positions,2) > 0_pInt + case ('memory_efficient') + memory_efficient = IO_intValue(line,positions,2) > 0_pInt !* Random seeding parameters case ('fixed_seed') @@ -322,7 +322,7 @@ subroutine numerics_init() write(6,'(a24,x,e8.1)') 'err_defgrad_tol: ',err_defgrad_tol write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,x,i8)') 'itmax: ',itmax - write(6,'(a24,x,L8)') 'fast_execution: ',fast_execution + write(6,'(a24,x,L8)') 'memory_efficient: ',memory_efficient write(6,*) !* Random seeding parameters @@ -377,10 +377,10 @@ subroutine numerics_init() if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(289) !* spectral parameters - if (err_div_tol <= 0.0_pReal) call IO_error(48) - if (err_defgrad_tol <= 0.0_pReal) call IO_error(48) - if (err_stress_tolrel <= 0.0_pReal) call IO_error(48) - if (itmax <= 1.0_pInt) call IO_error(48) + if (err_div_tol <= 0.0_pReal) call IO_error(49) + if (err_defgrad_tol <= 0.0_pReal) call IO_error(49) + if (err_stress_tolrel <= 0.0_pReal) call IO_error(49) + if (itmax <= 1.0_pInt) call IO_error(49) if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!' endsubroutine