223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

extended IO to cope with different name for solverJob and Model 

polishing, added error codes
added FFTW library files
This commit is contained in:
Martin Diehl 2011-02-21 14:37:38 +00:00
parent 34336bc659
commit 8dd1a694a3
10 changed files with 169 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
code/IO.f90
View File

@ -168,14 +168,14 @@ end function
if (FEsolver == 'Abaqus') then if (FEsolver == 'Abaqus') then
open(unit+1,status='old',err=100,& open(unit+1,status='old',err=100,&
file=trim(getSolverWorkingDirectoryName())//& file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension) trim(getModelName())//InputFileExtension)
open(unit,err=100,file=trim(getSolverWorkingDirectoryName())//& 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 IO_open_inputFile = IO_abaqus_assembleInputFile(unit,unit+1) ! strip comments and concatenate any "include"s
close(unit+1) close(unit+1)
else else
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//& open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension) trim(getModelName())//InputFileExtension)
IO_open_inputFile = .true. IO_open_inputFile = .true.
endif endif
@ -1142,10 +1142,12 @@ endfunction
case (45) case (45)
msg = 'error opening spectral loadcase' msg = 'error opening spectral loadcase'
case (46) case (46)
msg = 'missing parameter in spectral loadcase'
case (47)
msg = 'mask consistency violated in spectral loadcase' msg = 'mask consistency violated in spectral loadcase'
case (47)
msg = 'negative time increment in spectral loadcase'
case (48) case (48)
msg = 'Non-positive increment in spectral loadcase'
case (49)
msg = 'Non-positive relative parameter for spectral method' msg = 'Non-positive relative parameter for spectral method'
case (50) case (50)
msg = 'Error writing constitutive output description' msg = 'Error writing constitutive output description'
@ -1155,6 +1157,10 @@ endfunction
msg = 'Cannot open input file' msg = 'Cannot open input file'
case (102) case (102)
msg = 'argument count error (mesh and loadcase) for mpie_spectral' 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) case (105)
msg = 'Error reading from ODF file' msg = 'Error reading from ODF file'
case (110) case (110)
@ -1308,6 +1314,9 @@ endfunction
case (700) case (700)
msg = 'Singular matrix in stress iteration' msg = 'Singular matrix in stress iteration'
case (800)
msg = 'Matrix inversion error'
! Error messages related to parsing of Abaqus input file ! Error messages related to parsing of Abaqus input file
case (900) case (900)
msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)' msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)'

BIN
code/libfftw3.a Normal file
View File

Binary file not shown.

BIN
code/libfftw3_threads.a Normal file
View File

Binary file not shown.

8
code/makefile
View File

@ -6,11 +6,13 @@
# Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads" and "make", "make install" is not needed # 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. # as long as the two library files are copied to the source code directory.
cpspectral.out: mpie_spectral.o CPFEM.a cpspectral.exe: 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 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 mpie_spectral.o: mpie_spectral.f90 CPFEM.o
ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90 ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90
CPFEM.a: CPFEM.o CPFEM.a: CPFEM.o
ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o

18
code/mpie_cpfem_abaqus_exp.f
View File

@ -24,7 +24,9 @@ character(len=4), parameter :: InputFileExtension = '.inp'
CONTAINS CONTAINS
!--------------------
subroutine mpie_interface_init() subroutine mpie_interface_init()
!--------------------
write(6,*) write(6,*)
write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>' write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>'
write(6,*) '$Id$' write(6,*) '$Id$'
@ -32,7 +34,9 @@ subroutine mpie_interface_init()
return return
end subroutine end subroutine
!--------------------
function getSolverWorkingDirectoryName() function getSolverWorkingDirectoryName()
!--------------------
use prec use prec
implicit none implicit none
character(1024) getSolverWorkingDirectoryName character(1024) getSolverWorkingDirectoryName
@ -40,10 +44,21 @@ function getSolverWorkingDirectoryName()
getSolverWorkingDirectoryName='' getSolverWorkingDirectoryName=''
CALL VGETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR ) CALL VGETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR )
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function end function
!--------------------
function getModelName()
!--------------------
character(1024) getModelName
getModelName = getSolverJobName()
end function
!--------------------
function getSolverJobName() function getSolverJobName()
!--------------------
use prec use prec
implicit none implicit none
@ -52,7 +67,6 @@ function getSolverJobName()
getSolverJobName='' getSolverJobName=''
CALL VGETJOBNAME(getSolverJobName , LENJOBNAME ) CALL VGETJOBNAME(getSolverJobName , LENJOBNAME )
! write(6,*) 'getSolverJobName', getSolverJobName
end function end function
END MODULE END MODULE

16
code/mpie_cpfem_abaqus_std.f
View File

@ -24,7 +24,9 @@ character(len=4), parameter :: InputFileExtension = '.inp'
CONTAINS CONTAINS
!--------------------
subroutine mpie_interface_init() subroutine mpie_interface_init()
!--------------------
write(6,*) write(6,*)
write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>' write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>'
write(6,*) '$Id$' write(6,*) '$Id$'
@ -32,7 +34,9 @@ subroutine mpie_interface_init()
return return
end subroutine end subroutine
!--------------------
function getSolverWorkingDirectoryName() function getSolverWorkingDirectoryName()
!--------------------
use prec use prec
implicit none implicit none
character(1024) getSolverWorkingDirectoryName character(1024) getSolverWorkingDirectoryName
@ -43,7 +47,19 @@ function getSolverWorkingDirectoryName()
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName ! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function end function
!--------------------
function getModelName()
!--------------------
character(1024) getModelName
getModelName = getSolverJobName()
end function
!--------------------
function getSolverJobName() function getSolverJobName()
!--------------------
use prec use prec
implicit none implicit none

11
code/mpie_cpfem_marc.f90
View File

@ -71,6 +71,17 @@ function getSolverWorkingDirectoryName()
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName ! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function end function
function getModelName()
use prec
implicit none
character(1024) getModelName
getModelName = getSolverJobName()
end function
function getSolverJobName() function getSolverJobName()
use prec use prec
implicit none implicit none

137
code/mpie_spectral.f90
View File

@ -30,7 +30,7 @@ program mpie_spectral
use math use math
use CPFEM, only: CPFEM_general, CPFEM_initAll use CPFEM, only: CPFEM_general, CPFEM_initAll
use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol,& 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 homogenization, only: materialpoint_sizeResults, materialpoint_results
!$ use OMP_LIB ! the openMP function library !$ use OMP_LIB ! the openMP function library
@ -78,7 +78,7 @@ program mpie_spectral
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(3,3) :: xinormdyad
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat 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(pInt), dimension(3) :: k_s
integer*8, dimension(2) :: plan_fft integer*8, dimension(2) :: plan_fft
@ -109,8 +109,9 @@ program mpie_spectral
! Reading the loadcase file and assign variables ! Reading the loadcase file and assign variables
path = getLoadcaseName() path = getLoadcaseName()
print*,'Loadcase: ',trim(path) print '(a,/,a)', 'Loadcase: ',trim(path)
print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) 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) 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 N_n = N_n+1
end select end select
enddo ! count all identifiers to allocate memory and do sanity check 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 enddo
101 N_Loadcases = N_l 101 N_Loadcases = N_l
@ -150,7 +149,7 @@ program mpie_spectral
read(unit,'(a1024)',END = 200) line read(unit,'(a1024)',END = 200) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
i = i + 1 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 do j = 1,maxNchunksInput,2
select case (IO_lc(IO_stringValue(line,posInput,j))) select case (IO_lc(IO_stringValue(line,posInput,j)))
case('l','velocitygrad') case('l','velocitygrad')
@ -179,20 +178,22 @@ program mpie_spectral
200 close(unit) 200 close(unit)
do i = 1, N_Loadcases do i = 1, N_Loadcases
if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(47,i) ! bc_mask consistency if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(46,i) ! bc_mask consistency
print '(a,/,3(3(f12.6,x)/))','L' ,math_transpose3x3(bc_velocityGrad(:,:,i)) if (bc_timeIncrement(i) < 0.0_pReal) call IO_error(47,i) ! negative time increment forbidden
print '(a,/,3(3(f12.6,x)/))','bc_stress',math_transpose3x3(bc_stress(:,:,i)) if (bc_steps(i) < 1_pInt) call IO_error(48,i) ! non-positive increments requested
print '(a,/,3(3(l,x)/))', 'bc_mask for velocitygrad',transpose(bc_mask(:,:,1,i)) print '(a,/,3(3(f12.6,x)/))','L:' ,math_transpose3x3(bc_velocityGrad(:,:,i))
print '(a,/,3(3(l,x)/))', 'bc_mask for stress' ,transpose(bc_mask(:,:,2,i)) print '(a,/,3(3(f12.6,x)/))','bc_stress:',math_transpose3x3(bc_stress(:,:,i))
print *,'time',bc_timeIncrement(i) print '(a,/,3(3(l,x)/))', 'bc_mask for velocitygrad:',transpose(bc_mask(:,:,1,i))
print *,'incs',bc_steps(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 *, '' print *, ''
enddo enddo
!read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90 !read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90
path = getSolverJobName() path = getModelName()
print*,'JobName: ',trim(path) print '(a,a)', 'GeomName: ',trim(path)
if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = path) if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = trim(path)//InputFileExtension)
rewind(unit) rewind(unit)
do do
@ -233,11 +234,11 @@ program mpie_spectral
enddo enddo
100 close(unit) 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,/,i4,i4,i4)','resolution a b c:', resolution
print '(a,/,f6.1,f6.1,f6.1)','dimension x y z', geomdimension print '(a,/,f6.1,f6.1,f6.1)','dimension x y z:', geomdimension
print *,'homogenization',homog print '(a,i4)','homogenization: ',homog
allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal 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 allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
@ -261,10 +262,12 @@ program mpie_spectral
c066 = c066 * wgt c066 = c066 * wgt
c0 = math_mandel66to3333(c066) c0 = math_mandel66to3333(c066)
call math_invert(6, c066, s066,i, errmatinv) call math_invert(6, c066, s066,i, errmatinv)
if(errmatinv) call IO_error(800) !Matrix inversion error
s0 = math_mandel66to3333(s066) s0 = math_mandel66to3333(s066)
!calculation of calculate gamma_hat field in the case of fast execution (needs a lot of memory) if(memory_efficient) then ! allocate just single fourth order tensor
if(fast_execution) then 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 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) do k = 1, resolution(3)
k_s(3) = k-1 k_s(3) = k-1
@ -280,32 +283,31 @@ program mpie_spectral
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 if (any(xi /= 0.0_pReal)) then
do l = 1,3; do m = 1,3 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 enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad))
else else
xinormdyad = 0.0_pReal xinormdyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif 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 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)) *& gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(0.5*xinormdyad(m,p)+0.5*xinormdyad(p,m)) (xinormdyad(m,p)+xinormdyad(p,m))
enddo; enddo; enddo; enddo enddo; enddo; enddo; 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 endif
! calculate xi for the calculation of divergence in Fourier space (middle frequency) ! calculate xi for the calculation of divergence in Fourier space (central frequency)
xi_middle(3) = 0.0_pReal !2D case xi_central(3) = 0.0_pReal !2D case
if(resolution(3) > 1) xi_middle(3) = real(resolution(3)/2, pReal)/geomdimension(3) !3D case if(resolution(3) > 1) xi_central(3) = real(resolution(3)/2, pReal)/geomdimension(3) !3D case
xi_middle(2) = real(resolution(2)/2, pReal)/geomdimension(2) xi_central(2) = real(resolution(2)/2, pReal)/geomdimension(2)
xi_middle(1) = real(resolution(1)/2, pReal)/geomdimension(1) 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 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) ! 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_with_nthreads(mpieNumThreadsInt)
call dfftw_plan_many_dft_r2c(plan_fft(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& 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),& 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 ! write header of output file
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())&
//'_'//trim(getLoadcase())&
//'.spectralOut',form='UNFORMATTED') //'.spectralOut',form='UNFORMATTED')
write(538), 'load',trim(getLoadcaseName()) write(538), 'load',trim(getLoadcaseName())
write(538), 'workingdir',trim(getSolverWorkingDirectoryName()) write(538), 'workingdir',trim(getSolverWorkingDirectoryName())
@ -436,10 +437,10 @@ program mpie_spectral
enddo; enddo enddo; enddo
err_div = 2 * err_div_tol err_div = 2 * err_div_tol
err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) 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(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,/,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 stress: ',err_stress, ' Tol. = ', err_stress_tol
print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,' Tol. = ', err_defgrad_tol*0.8 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 if(err_stress < err_stress_tol*0.8) then
calcmode = 1_pInt calcmode = 1_pInt
endif endif
@ -457,10 +458,11 @@ program mpie_spectral
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt 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,:,:),& defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,& temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
CPFEM_mode = 2_pInt
workfft(i,j,k,:,:) = pstress workfft(i,j,k,:,:) = pstress
cstress_av = cstress_av + math_mandel6to33(cstress) cstress_av = cstress_av + math_mandel6to33(cstress)
enddo; enddo; enddo 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))) sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:) + (workfft(2,1,1,m,:))*img)))
enddo 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) 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 err_div = err_div/sigma0 !weighting of error
if(fast_execution) then ! fast execution with stored gamma_hat if(memory_efficient) then ! memory saving version, in-time calculation of 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
do k = 1, resolution(3) do k = 1, resolution(3)
k_s(3) = k-1 k_s(3) = k-1
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) 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 xi(3) = 0.0_pReal !for the 2D case
if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) !3D 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(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 if (any(xi(:) /= 0.0_pReal)) then
do l = 1,3; do m = 1,3 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 enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad))
else else
xinormdyad = 0.0_pReal xinormdyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif 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 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))*& gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*&
(0.5*xinormdyad(m,p) +0.5*xinormdyad(p,m)) (xinormdyad(m,p) +xinormdyad(p,m))
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
do m = 1,3; do n = 1,3 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,:,:)& 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-1,j,k,:,:) = real (temp33_Complex)
workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex)
enddo; enddo; enddo 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 endif
workfft(1,1,1,:,:) = defgrad_av - math_I3 !zero frequency (real part) 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_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel !accecpt relativ error specified
err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) 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 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 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 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. 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 calcmode = 0_pInt
@ -552,10 +553,10 @@ program mpie_spectral
write(538) materialpoint_results(:,1,:) !write to output file write(538) materialpoint_results(:,1,:) !write to output file
print '(a,x,f12.7)' , ' Determinant of Deformation Aim:', math_det3x3(defgradAim) 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 Aim: ',math_transpose3x3(defgradAim)
print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av) 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,/,3(3(f10.4,x)/))', ' Cauchy Stress / MPa: ',math_transpose3x3(cstress_av)/1.e6
print '(A)', '************************************************************' print '(A)', '************************************************************'
enddo ! end looping over steps in current loadcase enddo ! end looping over steps in current loadcase
enddo ! end looping over loadcases enddo ! end looping over loadcases

31
code/mpie_spectral_interface.f90
View File

@ -59,28 +59,43 @@ function getSolverJobName()
implicit none 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) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \
integer(pInt) posExt,posSep integer(pInt) posExt,posSep
getSolverJobName = '' getModelName = ''
call getarg(1,outName) call getarg(1,outName)
posExt = scan(outName,'.',back=.true.) posExt = scan(outName,'.',back=.true.)
posSep = scan(outName,pathSep,back=.true.) posSep = scan(outName,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(outName)+1 ! no extension present 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) call getcwd(cwd)
getSolverJobName = rectifyPath(trim(cwd)//'/'//getSolverJobName) getModelName = rectifyPath(trim(cwd)//'/'//getModelName)
else else
getSolverJobName = rectifyPath(getSolverJobName) getModelName = rectifyPath(getModelName)
endif endif
getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),& getModelName = makeRelativePath(getSolverWorkingDirectoryName(),&
getSolverJobName) getModelName)
return return
endfunction endfunction

30
code/numerics.f90
View File

@ -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_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_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 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 integer(pInt) itmax , & ! maximum number of iterations
@ -94,7 +94,7 @@ subroutine numerics_init()
character(len=1024) line character(len=1024) line
! OpenMP variable ! OpenMP variable
!$ character(len=2) mpieNumThreadsString !enironment variable MPIE_NUMTHREADS !$ character(len=4) mpieNumThreadsString !environment variable MPIE_NUMTHREADS
write(6,*) write(6,*)
write(6,*) '<<<+- numerics init -+>>>' write(6,*) '<<<+- numerics init -+>>>'
@ -142,11 +142,11 @@ subroutine numerics_init()
volDiscrPow_RGC = 5.0 volDiscrPow_RGC = 5.0
!* spectral parameters: !* spectral parameters:
err_div_tol = 1.0e-4 err_div_tol = 1.0e-4 ! proposed by Suquet, less strict criteria are usefull, e.g. 5e-3
err_defgrad_tol = 1.0e-3 err_defgrad_tol = 1.0e-3 ! relative tolerance for fullfillment of average deformation gradient (is usually passively fullfilled)
err_stress_tolrel = 0.01 err_stress_tolrel = 0.01 ! relative tolerance for fullfillment of stress BC
itmax = 20_pInt itmax = 20_pInt ! Maximum iteration number
fast_execution = .false. memory_efficient = .true. ! Precalculate Gamma-operator (81 double per point)
!* Random seeding parameters: added <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
fixedSeed = 0_pInt fixedSeed = 0_pInt
@ -154,7 +154,7 @@ subroutine numerics_init()
!* determin number of threads from environment variable MPIE_NUM_THREADS !* determin number of threads from environment variable MPIE_NUM_THREADS
!$ call GetEnv('MPIE_NUM_THREADS',mpieNumThreadsString) ! get 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... !$ 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 !$ 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) err_stress_tolrel = IO_floatValue(line,positions,2)
case ('itmax') case ('itmax')
itmax = IO_intValue(line,positions,2) itmax = IO_intValue(line,positions,2)
case ('fast_execution') case ('memory_efficient')
fast_execution = IO_intValue(line,positions,2) > 0_pInt memory_efficient = IO_intValue(line,positions,2) > 0_pInt
!* Random seeding parameters !* Random seeding parameters
case ('fixed_seed') 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_defgrad_tol: ',err_defgrad_tol
write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,x,i8)') 'itmax: ',itmax 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,*) write(6,*)
!* Random seeding parameters !* Random seeding parameters
@ -377,10 +377,10 @@ subroutine numerics_init()
if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(289) if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(289)
!* spectral parameters !* spectral parameters
if (err_div_tol <= 0.0_pReal) 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(48) if (err_defgrad_tol <= 0.0_pReal) call IO_error(49)
if (err_stress_tolrel <= 0.0_pReal) call IO_error(48) if (err_stress_tolrel <= 0.0_pReal) call IO_error(49)
if (itmax <= 1.0_pInt) call IO_error(48) if (itmax <= 1.0_pInt) call IO_error(49)
if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!' if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!'
endsubroutine endsubroutine