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

restructured algorithm, initialization now not longer within increments, lot of small improvements/polishing 

makefile now calls compiler with lot of warning flags
This commit is contained in:
Martin Diehl 2012-02-10 11:59:59 +00:00
parent 37ac7bf1b4
commit 1cc2315954
3 changed files with 516 additions and 387 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

582
code/DAMASK_spectral.f90
View File

@ -49,7 +49,7 @@ program DAMASK_spectral
use math use math
use kdtree2_module use kdtree2_module
use CPFEM, only: CPFEM_general, CPFEM_initAll use CPFEM, only: CPFEM_general, CPFEM_initAll
use FEsolving, only: restartWrite, restartReadInc use FEsolving, only: restartWrite, restartInc
use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, & use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, &
memory_efficient, update_gamma, & memory_efficient, update_gamma, &
simplified_algorithm, divergence_correction, & simplified_algorithm, divergence_correction, &
@ -115,13 +115,23 @@ program DAMASK_spectral
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: pstress, pstress_av, & real(pReal), dimension(3,3) :: pstress, pstress_av, &
defgradAim = math_I3, defgradAimOld = math_I3,& defgradAim = math_I3, defgradAimOld = math_I3,&
mask_stress, mask_defgrad, fDot, & mask_stress, mask_defgrad, deltaF, &
pstress_av_lab, defgradAim_lab, defgrad_av_lab ! quantities rotated to other coordinate system pstress_av_lab, defgradAim_lab, defgrad_av_lab ! quantities rotated to other coordinate system
real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current = 0.0_pReal, s_prev, c_prev ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current = 0.0_pReal, s_prev, c_prev,& ! stiffness and compliance
s0_reference
real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6) :: cstress ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde, c0_66, s0_66 ! small strain stiffness
real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation real(pReal), dimension(9,9) :: s_prev99, c_prev99, c0_99, s0_99 ! compliance and stiffness in matrix notation
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
real(pReal), dimension(6,6) :: mask_inversion = reshape([&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal],&
[ 6_pInt, 6_pInt])
real(pReal), dimension(3,3,3,3) :: temp_3333 = 0.0_pReal
integer(pInt) :: size_reduced = 0.0_pReal ! number of stress BCs integer(pInt) :: size_reduced = 0.0_pReal ! number of stress BCs
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -129,8 +139,8 @@ program DAMASK_spectral
type(C_PTR) :: tensorField, tau ! fields in real an fourier space type(C_PTR) :: tensorField, tau ! fields in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: tensorField_real ! fields in real space (pointer) real(pReal), dimension(:,:,:,:,:), pointer :: tensorField_real ! fields in real space (pointer)
real(pReal), dimension(:,:,:,:,:), pointer :: tau_real real(pReal), dimension(:,:,:,:,:), pointer :: tau_real
complex(pReal), dimension(:,:,:,:,:), pointer :: tensorField_complex ! fields in fourier space (pointer) complex(pReal), dimension(:,:,:,:,:), pointer :: tensorField_fourier ! fields in fourier space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: tau_complex complex(pReal), dimension(:,:,:,:,:), pointer :: tau_fourier
real(pReal), dimension(:,:,:,:,:), allocatable :: defgrad, defgradold real(pReal), dimension(:,:,:,:,:), allocatable :: defgrad, defgradold
real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature real(pReal), dimension(:,:,:), allocatable :: temperature
@ -145,7 +155,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc ! elapsed time, begin of interval, time interval real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval
real(pReal) :: guessmode, err_div, err_stress, err_stress_tol real(pReal) :: guessmode, err_div, err_stress, err_stress_tol
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3) :: temp3_Complex complex(pReal), dimension(3) :: temp3_Complex
@ -153,8 +163,8 @@ program DAMASK_spectral
real(pReal), dimension(3,3) :: temp33_Real real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, m, n, p, errorID integer(pInt) :: i, j, k, l, m, n, p, errorID
integer(pInt) :: N_Loadcases, loadcase, inc, iter, ielem, CPFEM_mode, & integer(pInt) :: N_Loadcases, loadcase, inc, iter, ielem, CPFEM_mode, &
ierr, notConvergedCounter = 0_pInt, totalIncsCounter = 0_pInt,& ierr, totalIncsCounter = 0_pInt,&
writtenRestart = 0_pInt notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
logical :: errmatinv logical :: errmatinv
real(pReal) :: defgradDet, correctionFactor real(pReal) :: defgradDet, correctionFactor
@ -164,25 +174,24 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!variables for additional output due to general debugging !variables for additional output due to general debugging
real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew, max_diag, max_offdiag
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations ! variables for additional output of divergence calculations
type(C_PTR) :: divergence, plan_divergence type(C_PTR) :: divergence, plan_divergence
real(pReal), dimension(:,:,:,:), pointer :: divergence_real real(pReal), dimension(:,:,:,:), pointer :: divergence_real
complex(pReal), dimension(:,:,:,:), pointer :: divergence_complex complex(pReal), dimension(:,:,:,:), pointer :: divergence_fourier
real(pReal), dimension(:,:,:,:), allocatable :: divergence_postProc real(pReal), dimension(:,:,:,:), allocatable :: divergence_postProc
real(pReal) :: p_hat_avg, p_real_avg,& real(pReal) :: pstress_av_L2, err_div_RMS, err_real_div_RMS,&
err_div_RMS, err_real_div_RMS,&
err_div_max, err_real_div_max,& err_div_max, err_real_div_max,&
max_div_error max_div_error
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables for debugging fft using a scalar field ! variables for debugging fft using a scalar field
type(C_PTR) :: scalarField_realPointer, scalarField_complexPointer,& type(C_PTR) :: scalarField_realC, scalarField_fourierC,&
plan_scalarField_forth, plan_scalarField_back plan_scalarField_forth, plan_scalarField_back
real(pReal), dimension(:,:,:), pointer :: scalarField_real complex(pReal), dimension(:,:,:), pointer :: scalarField_real
complex(pReal), dimension(:,:,:), pointer :: scalarField_complex complex(pReal), dimension(:,:,:), pointer :: scalarField_fourier
integer(pInt) :: row, column integer(pInt) :: row, column
!################################################################################################## !##################################################################################################
@ -250,7 +259,7 @@ program DAMASK_spectral
do k = 1_pInt,9_pInt do k = 1_pInt,9_pInt
if (temp_maskVector(k)) temp_valueVector(k) = IO_floatValue(line,positions,j+k) if (temp_maskVector(k)) temp_valueVector(k) = IO_floatValue(line,positions,j+k)
enddo enddo
bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,(/3,3/))) bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,[ 3,3]))
bc(loadcase)%deformation = math_plain9to33(temp_valueVector) bc(loadcase)%deformation = math_plain9to33(temp_valueVector)
case('p','pk1','piolakirchhoff','stress') case('p','pk1','piolakirchhoff','stress')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
@ -260,7 +269,7 @@ program DAMASK_spectral
if (bc(loadcase)%maskStressVector(k)) temp_valueVector(k) =& if (bc(loadcase)%maskStressVector(k)) temp_valueVector(k) =&
IO_floatValue(line,positions,j+k) ! assign values for the bc(loadcase)%stress matrix IO_floatValue(line,positions,j+k) ! assign values for the bc(loadcase)%stress matrix
enddo enddo
bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,(/3,3/))) bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,[ 3,3]))
bc(loadcase)%stress = math_plain9to33(temp_valueVector) bc(loadcase)%stress = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt) bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt)
@ -359,7 +368,7 @@ program DAMASK_spectral
! sanity checks of geometry parameters ! sanity checks of geometry parameters
if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))& if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))&
call IO_error(error_ID = 45_pInt) call IO_error(error_ID = 45_pInt)
if (any(geomdim<=0.0_pReal)) stop if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 102_pInt)
if(mod(res(1),2_pInt)/=0_pInt .or.& if(mod(res(1),2_pInt)/=0_pInt .or.&
mod(res(2),2_pInt)/=0_pInt .or.& mod(res(2),2_pInt)/=0_pInt .or.&
(mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt))& (mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt))&
@ -411,9 +420,9 @@ program DAMASK_spectral
print '(a)','deformation gradient rate:' print '(a)','deformation gradient rate:'
endif endif
write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),transpose(bc(loadcase)%maskDeformation)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation))
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',&
1e-9*merge(math_transpose33(bc(loadcase)%stress),reshape(spread(DAMASK_NaN,1,9),(/3,3/))& 1e-9*merge(math_transpose33(bc(loadcase)%stress),reshape(spread(DAMASK_NaN,1,9),[ 3,3])&
,transpose(bc(loadcase)%maskStress)) ,transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',&
@ -426,10 +435,10 @@ program DAMASK_spectral
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only
if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. & if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/)))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 38_pInt ! no rotation is allowed by stress BC errorID = 38_pInt ! no rotation is allowed by stress BC
if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose33(bc(loadcase)%rotation))& if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose33(bc(loadcase)%rotation))&
-math_I3) > reshape(spread(rotation_tol,1,9),(/3,3/)))& -math_I3) > reshape(spread(rotation_tol,1,9),[ 3,3]))&
.or. abs(math_det33(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)& .or. abs(math_det33(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)&
errorID = 46_pInt ! given rotation matrix contains strain errorID = 46_pInt ! given rotation matrix contains strain
if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment
@ -446,56 +455,9 @@ program DAMASK_spectral
debugFFTW = iand(debug_spectral,debug_spectralFFTW) > 0_pInt debugFFTW = iand(debug_spectral,debug_spectralFFTW) > 0_pInt
!################################################################################################## !##################################################################################################
! Loop over loadcases defined in the loadcase file ! initialization
!################################################################################################## !##################################################################################################
do loadcase = 1_pInt, N_Loadcases
time0 = time ! loadcase start time
if (bc(loadcase)%followFormerTrajectory) then ! continue to guess along former trajectory where applicable
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif
!--------------------------------------------------------------------------------------------------
! arrays for mixed boundary conditions
mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation)
mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress)
size_reduced = count(bc(loadcase)%maskStressVector)
allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal
allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal
timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
fDot = bc(loadcase)%deformation ! only valid for given fDot. will be overwritten later in case L is given
!##################################################################################################
! loop oper incs defined in input file for current loadcase
!##################################################################################################
do inc = 1_pInt, bc(loadcase)%incs
!--------------------------------------------------------------------------------------------------
! forwarding time
if (bc(loadcase)%logscale == 1_pInt) then ! logarithmic scale
if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal))
endif
else ! not-1st loadcase of logarithmic scale
timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**(real( inc,pReal)/&
real(bc(loadcase)%incs ,pReal))&
-(1.0_pReal + bc(loadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(bc(loadcase)%incs ,pReal)) )
endif
endif
time = time + timeinc
totalIncsCounter = totalIncsCounter + 1_pInt
!##################################################################################################
! initialization start after forwarding to restart step
!##################################################################################################
if(totalIncsCounter == restartReadInc+1_pInt) then ! Initialize values
guessmode = 0.0_pReal ! no old values
allocate (defgrad ( res(1), res(2),res(3),3,3)); defgrad = 0.0_pReal allocate (defgrad ( res(1), res(2),res(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold ( res(1), res(2),res(3),3,3)); defgradold = 0.0_pReal allocate (defgradold ( res(1), res(2),res(3),3,3)); defgradold = 0.0_pReal
allocate (coordinates( res(1), res(2),res(3),3)); coordinates = 0.0_pReal allocate (coordinates( res(1), res(2),res(3),3)); coordinates = 0.0_pReal
@ -503,83 +465,32 @@ program DAMASK_spectral
allocate (xi (3,res1_red,res(2),res(3))); xi = 0.0_pReal allocate (xi (3,res1_red,res(2),res(3))); xi = 0.0_pReal
tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer(tensorField, tensorField_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation call c_f_pointer(tensorField, tensorField_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation
call c_f_pointer(tensorField, tensorField_complex, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation call c_f_pointer(tensorField, tensorField_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details) ! init fields to no deformation
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102) ! check for correct precision in C ielem = 0_pInt
#ifdef _OPENMP
if(DAMASK_NumThreadsInt > 0_pInt) then
ierr = fftw_init_threads()
if (ierr == 0_pInt) call IO_error(error_ID = 104_pInt)
call fftw_plan_with_nthreads(DAMASK_NumThreadsInt)
endif
#endif
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
!--------------------------------------------------------------------------------------------------
! creating plans
plan_stress = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),9,& ! dimensions , length in each dimension in reversed order
tensorField_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
tensorField_complex,(/res(3),res(2) ,res1_red/),&
1, res(3)*res(2)* res1_red,fftw_planner_flag)
plan_correction =fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),9,&
tensorField_complex,(/res(3),res(2) ,res1_red/),&
1, res(3)*res(2)* res1_red,&
tensorField_real,(/res(3),res(2) ,res(1)+2_pInt/),&
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
!--------------------------------------------------------------------------------------------------
! depending on (debug) options, allocate more memory and create additional plans
if (.not. simplified_algorithm) then
print*, 'using polarization field based algorithm'
tau = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T))
call c_f_pointer(tau, tau_real, [ res(1)+2_pInt,res(2),res(3),3,3])
call c_f_pointer(tau, tau_complex, [ res1_red, res(2),res(3),3,3])
plan_tau = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),9,&
tau_real,(/res(3),res(2) ,res(1)+2_pInt/),&
1, res(3)*res(2)*(res(1)+2_pInt),&
tau_complex,(/res(3),res(2) ,res1_red/),&
1, res(3)*res(2)* res1_red,fftw_planner_flag)
endif
if (debugDivergence) then
divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T))
call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3])
call c_f_pointer(divergence, divergence_complex, [ res1_red, res(2),res(3),3])
allocate (divergence_postProc(res(1),res(2),res(3),3)); divergence_postProc= 0.0_pReal
plan_divergence = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),3,&
divergence_complex,(/res(3),res(2) ,res1_red/),&
1, res(3)*res(2)* res1_red,&
divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),&
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
endif
if (debugFFTW) then
scalarField_realPointer = fftw_alloc_complex(int(res(1) *res(2)*res(3),C_SIZE_T)) ! do not do an inplace transform
scalarField_complexPointer = fftw_alloc_complex(int(res1_red*res(2)*res(3),C_SIZE_T))
call c_f_pointer(scalarField_realPointer, scalarField_real, [res(1), res(2),res(3)])
call c_f_pointer(scalarField_complexPointer, scalarField_complex, [res1_red,res(2),res(3)])
plan_scalarField_forth = fftw_plan_dft_r2c_3d(res(3),res(2),res(1),& !reversed order
scalarField_real,scalarField_complex,fftw_planner_flag)
plan_scalarField_back = fftw_plan_dft_c2r_3d(res(3),res(2),res(1),& !reversed order
scalarField_complex,scalarField_real,fftw_planner_flag)
endif
if (debugGeneral) print '(a)' , 'FFTW initialized'
!--------------------------------------------------------------------------------------------------
! calculate initial deformation
if (restartReadInc==0_pInt) then ! not restarting, no deformation at the beginning
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
defgrad(i,j,k,1:3,1:3) = math_I3 defgrad(i,j,k,1:3,1:3) = math_I3
defgradold(i,j,k,1:3,1:3) = math_I3 defgradold(i,j,k,1:3,1:3) = math_I3
coordinates(i,j,k,1:3) = geomdim/real(res, pReal)*[i,j,k] - geomdim/real(2_pInt*res,pReal)
call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c_current = c_current + dPdF
enddo; enddo; enddo enddo; enddo; enddo
else ! using old values from file
c0_reference = c_current * wgt ! linear reference material stiffness
c0_66 = math_Mandel3333to66(c0_reference)
call math_invert(6_pInt, c0_66, s0_66, i, errmatinv) ! invert in mandel notation
if(errmatinv) call IO_error(error_ID=800_pInt)
s0_reference = math_Mandel66to3333(s0_66)
!--------------------------------------------------------------------------------------------------
! possible restore deformation gradient from saved state
if (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',& if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',&
restartReadInc ,' from file' restartInc - 1_pInt,' from file'
if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(defgrad))) then trim(getSolverJobName()),size(defgrad))) then
read (777,rec=1) defgrad read (777,rec=1) defgrad
@ -592,22 +503,10 @@ program DAMASK_spectral
enddo; enddo; enddo enddo; enddo; enddo
defgradAim = defgradAim * wgt defgradAim = defgradAim * wgt
defgradAimOld = defgradAim defgradAimOld = defgradAim
guessmode=0.0_pInt
endif endif
call deformed_fft(res,geomdim,defgradAimOld,1.0_pReal,defgrad,coordinates) ! calculate current coordinates
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c_current = c_current + dPdF
enddo; enddo; enddo
c0_reference = c_current * wgt ! linear reference material stiffness
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies
if (debugGeneral) print '(a)' , 'first call to CPFEM_general finished'
do k = 1_pInt, res(3) do k = 1_pInt, res(3)
k_s(3) = k - 1_pInt k_s(3) = k - 1_pInt
if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3)
@ -619,12 +518,6 @@ program DAMASK_spectral
xi(1:3,i,j,k) = real(k_s, pReal)/geomdim xi(1:3,i,j,k) = real(k_s, pReal)/geomdim
enddo; enddo; enddo enddo; enddo; enddo
xi(1,res1_red-cutting_freq(1):res1_red , 1:res(2) , 1:res(3)) = 0.0_pReal
xi(2,1:res1_red, res(2)/2_pInt+1_pInt-cutting_freq(2):res(2)/2_pInt+1_pInt+cutting_freq(2),&
1:res(3)) = 0.0_pReal
xi(3,1:res1_red, 1:res(2) ,&
res(3)/2_pInt+1_pInt-cutting_freq(3):res(3)/2_pInt+1_pInt+cutting_freq(3)) = 0.0_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate the gamma operator ! calculate the gamma operator
if(memory_efficient) then ! allocate just single fourth order tensor if(memory_efficient) then ! allocate just single fourth order tensor
@ -633,34 +526,98 @@ program DAMASK_spectral
allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then
do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt do m = 1_pInt ,3_pInt; do p = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) xiDyad(m,p) = xi(m, i,j,k)*xi(p, i,j,k)
enddo; enddo enddo; enddo
temp33_Real = math_inv33(math_mul3333xx33(c0_reference, xiDyad)) ! do l = 1_pInt,3_pInt
! do n = 1_pInt,3_pInt
! temp33_Real(l,n) = sum(c0_reference(l,1:3,n,1:3)*xiDyad(1:3,1:3))
! enddo; enddo
temp33_Real= math_mul3333xx33(c0_reference,xiDyad)
temp33_Real = math_inv33(temp33_Real)
else else
xiDyad = 0.0_pReal xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
endif endif
! if (k==res(3)/2 .or. k==res(3)/2+2 .or.&
! j==res(2)/2 .or. j==res(2)/2+2 .or.&
! i==res(1)/2 .or. i==res(1)/2+2) &
! gamma_hat(1,1,1,1:3,1:3,1:3,1:3) = s0_reference
do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt
gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *& gamma_hat(i,j,k, l,m,n,p) = 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(xiDyad(m,p)+xiDyad(p,m)) (xiDyad(m,p)+xiDyad(p,m))
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
enddo; enddo; enddo enddo; enddo; enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! empirical factor for making divergence resolution and dimension indpendent ! general initialization of fftw (see manual on fftw.org for more details)
divergence_correction =.false. if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=108_pInt) ! check for correct precision in C
if (divergence_correction) then #ifdef _OPENMP
if (res(3) == 1_pInt) then if(DAMASK_NumThreadsInt > 0_pInt) then
correctionFactor = minval(geomdim(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct? ierr = fftw_init_threads()
else if (ierr == 0_pInt) call IO_error(error_ID = 109_pInt)
correctionFactor = minval(geomdim(1:3))*wgt**(-1.0_pReal/4.0_pReal) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency call fftw_plan_with_nthreads(DAMASK_NumThreadsInt)
endif endif
else #endif
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
!--------------------------------------------------------------------------------------------------
! creating plans
plan_stress = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order
tensorField_real,[ res(3),res(2) ,res(1)+2_pInt],& ! input data , physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
tensorField_fourier,[ res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,fftw_planner_flag)
plan_correction =fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],9,&
tensorField_fourier,[ res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,&
tensorField_real,[ res(3),res(2) ,res(1)+2_pInt],&
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
!--------------------------------------------------------------------------------------------------
! depending on (debug) options, allocate more memory and create additional plans
if (.not. simplified_algorithm) then
print*, 'using polarization field based algorithm'
tau = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T))
call c_f_pointer(tau, tau_real, [ res(1)+2_pInt,res(2),res(3),3,3])
call c_f_pointer(tau, tau_fourier, [ res1_red, res(2),res(3),3,3])
plan_tau = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,&
tau_real,[ res(3),res(2) ,res(1)+2_pInt],&
1, res(3)*res(2)*(res(1)+2_pInt),&
tau_fourier,[ res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,fftw_planner_flag)
endif
if (debugDivergence) then
divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T))
call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3])
call c_f_pointer(divergence, divergence_fourier, [ res1_red, res(2),res(3),3])
allocate (divergence_postProc(res(1),res(2),res(3),3)); divergence_postProc= 0.0_pReal
plan_divergence = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],3,&
divergence_fourier,[ res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,&
divergence_real,[ res(3),res(2) ,res(1)+2_pInt],&
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
endif
if (debugFFTW) then
scalarField_realC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! do not do an inplace transform
scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T))
call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)])
call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)])
plan_scalarField_forth = fftw_plan_dft_3d(res(3),res(2),res(1),& !reversed order
scalarField_real,scalarField_fourier,-1,fftw_planner_flag)
plan_scalarField_back = fftw_plan_dft_3d(res(3),res(2),res(1),& !reversed order
scalarField_fourier,scalarField_real,+1,fftw_planner_flag)
endif
if (debugGeneral) print '(a)' , 'FFTW initialized'
!--------------------------------------------------------------------------------------------------
! do not correct divergence criterion (usefull to kill dimension and resolution dependenc)
correctionFactor = 1.0_pReal correctionFactor = 1.0_pReal
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
@ -676,55 +633,87 @@ program DAMASK_spectral
write(538) 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase write(538) 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538) 'times', bc(1:N_Loadcases)%time ! one entry per loadcase write(538) 'times', bc(1:N_Loadcases)%time ! one entry per loadcase
write(538) 'logscales', bc(1:N_Loadcases)%logscale write(538) 'logscales', bc(1:N_Loadcases)%logscale
bc(1)%incs = bc(1)%incs + 1_pInt ! additional for zero deformation
write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase
bc(1)%incs = bc(1)%incs - 1_pInt write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc
write(538) 'startingIncrement', restartReadInc ! start with writing out the previous inc
write(538) 'eoh' ! end of header write(538) 'eoh' ! end of header
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
if (debugGeneral) print '(a)' , 'Header of result file written out' if (debugGeneral) print '(a)' , 'Header of result file written out'
!##################################################################################################
! Loop over loadcases defined in the loadcase file
!##################################################################################################
do loadcase = 1_pInt, N_Loadcases
time0 = time ! loadcase start time
if (bc(loadcase)%followFormerTrajectory .and. &
(restartInc < totalIncsCounter .or. &
restartInc > totalIncsCounter+bc(loadcase)%incs) ) then ! continue to guess along former trajectory where applicable
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif endif
if(totalIncsCounter > restartReadInc) then ! Do calculations (otherwise just forwarding) !--------------------------------------------------------------------------------------------------
if(bc(loadcase)%restartFrequency>0_pInt) & ! arrays for mixed boundary conditions
restartWrite = ( mod(inc - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation)
if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress)
fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim) size_reduced = count(bc(loadcase)%maskStressVector)
allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal
allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal
!##################################################################################################
! loop oper incs defined in input file for current loadcase
!##################################################################################################
do inc = 1_pInt, bc(loadcase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
!--------------------------------------------------------------------------------------------------
! forwarding time
timeinc_old = timeinc
if (bc(loadcase)%logscale == 0_pInt) then ! linear scale
timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else
if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st loadcase of logarithmic scale
timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal))
endif
else ! not-1st loadcase of logarithmic scale
timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**(real( inc,pReal)/&
real(bc(loadcase)%incs ,pReal))&
-(1.0_pReal + bc(loadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(bc(loadcase)%incs ,pReal)) )
endif
endif
time = time + timeinc
if (bc(loadcase)%velGradApplied) then ! calculate deltaF from given L and current F
deltaF = timeinc * mask_defgrad * math_mul33x33(bc(loadcase)%deformation, defgradAim)
else ! deltaF = fDot *timeinc where applicable
deltaF = timeinc * mask_defgrad * bc(loadcase)%deformation
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system ! winding forward of deformation aim in loadcase system
temp33_Real = defgradAim temp33_Real = defgradAim
defgradAim = defgradAim & defgradAim = defgradAim &
+ guessmode * mask_stress * (defgradAim - defgradAimOld) & + guessmode * mask_stress * (defgradAim - defgradAimOld)*timeinc/timeinc_old &
+ mask_defgrad * fDot * timeinc + deltaF
defgradAimOld = temp33_Real defgradAimOld = temp33_Real
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update local deformation gradient ! update local deformation gradient
if (any(bc(loadcase)%rotation/=math_I3)) then ! lab and loadcase coordinate system are NOT the same deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = defgrad(i,j,k,1:3,1:3) temp33_Real = defgrad(i,j,k,1:3,1:3)
if (bc(loadcase)%velGradApplied) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33(bc(loadcase)%deformation,&
math_rotate_forward33(defgradold(i,j,k,1:3,1:3),bc(loadcase)%rotation))
defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing... + guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing...
+ math_rotate_backward33((1.0_pReal-guessmode) * mask_defgrad * fDot,& *timeinc/timeinc_old &
bc(loadcase)%rotation) *timeinc ! apply the prescribed value where deformation is given if not guessing + (1.0_pReal-guessmode) * deltaF ! if not guessing, use prescribed average deformation where applicable
defgradold(i,j,k,1:3,1:3) = temp33_Real defgradold(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
else ! one coordinate system for lab and loadcase, save some multiplications
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = defgrad(i,j,k,1:3,1:3)
if (bc(loadcase)%velGradApplied) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33(bc(loadcase)%deformation,defgradold(i,j,k,1:3,1:3))
defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing...
+ (1.0_pReal-guessmode) * mask_defgrad * fDot * timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate reduced compliance ! calculate reduced compliance
@ -742,7 +731,7 @@ program DAMASK_spectral
c_reduced(k,j) = c_prev99(n,m) c_reduced(k,j) = c_prev99(n,m)
endif; enddo; endif; enddo endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=799) if(errmatinv) call IO_error(error_ID=800_pInt)
s_prev99 = 0.0_pReal ! build full compliance s_prev99 = 0.0_pReal ! build full compliance
k = 0_pInt k = 0_pInt
do n = 1_pInt,9_pInt do n = 1_pInt,9_pInt
@ -761,14 +750,6 @@ program DAMASK_spectral
! report begin of new increment ! report begin of new increment
print '(a)', '##################################################################' print '(a)', '##################################################################'
print '(A,I5.5,A,es12.6)', 'Increment ', totalIncsCounter, ' Time ',time print '(A,I5.5,A,es12.6)', 'Increment ', totalIncsCounter, ' Time ',time
if (restartWrite ) then
print '(A)', 'writing converged results of previous increment for restart'
if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! writing deformation gradient field to file
write (777,rec=1) defgrad
close (777)
endif
writtenRestart=totalIncsCounter-1_pInt
endif
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
@ -827,7 +808,6 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! copy one component of the stress field to to a single FT and check for mismatch ! copy one component of the stress field to to a single FT and check for mismatch
if (debugFFTW) then if (debugFFTW) then
scalarField_real = 0.0_pReal
row = (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 row = (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1
column = (mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt column = (mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt
scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component
@ -843,7 +823,7 @@ program DAMASK_spectral
= tensorField_real(i,j,k,1:3,1:3) & = tensorField_real(i,j,k,1:3,1:3) &
- math_mul3333xx33(c0_reference,defgrad(i,j,k,1:3,1:3)) - math_mul3333xx33(c0_reference,defgrad(i,j,k,1:3,1:3))
enddo; enddo; enddo enddo; enddo; enddo
call fftw_execute_dft_r2c(plan_tau,tau_real,tau_complex) call fftw_execute_dft_r2c(plan_tau,tau_real,tau_fourier)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -854,12 +834,37 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing the FT because it simplifies calculation of average stress in real space also ! doing the FT because it simplifies calculation of average stress in real space also
call fftw_execute_dft_r2c(plan_stress,tensorField_real,tensorField_complex) call fftw_execute_dft_r2c(plan_stress,tensorField_real,tensorField_fourier)
pstress_av_lab = real(tensorField_complex(1,1,1,1:3,1:3),pReal)*wgt
pstress_av_lab = real(tensorField_fourier(1,1,1,1:3,1:3),pReal)*wgt
pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation) pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation)
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa:',&
math_transpose33(pstress_av)/1.e6 math_transpose33(pstress_av)/1.e6
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 FT results
if (debugFFTW) then
call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier)
print '(a,i1,1x,i1)', 'checking FT results of compontent ', row, column
print '(a,2(es10.4,1x))', 'max FT relative error ',&
maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
tensorField_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), &
maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
tensorField_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3))))
endif
!--------------------------------------------------------------------------------------------------
! removing highest frequencies
tensorField_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal)
tensorField_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)>1_pInt) &
tensorField_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
@ -884,61 +889,50 @@ program DAMASK_spectral
print '(a)', '' print '(a)', ''
print '(a)', '... calculating equilibrium with spectral method .................' print '(a)', '... calculating equilibrium with spectral method .................'
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 FT results
if (debugFFTW) then
call fftw_execute_dft_r2c(plan_scalarField_forth,scalarField_real,scalarField_complex)
print '(a,i1,1x,i1)', 'checking FT results of compontent ', row, column
print '(a,2(es10.4,1x))', 'max FT relative error ',&
maxval( real((scalarField_complex(1:res1_red,1:res(2),1:res(3))-&
tensorField_complex(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_complex(1:res1_red,1:res(2),1:res(3)))), &
maxval(aimag((scalarField_complex(1:res1_red,1:res(2),1:res(3))-&
tensorField_complex(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_complex(1:res1_red,1:res(2),1:res(3))))
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space ! calculating RMS divergence criterion in Fourier space
p_hat_avg = sqrt(maxval (math_eigenvalues33(math_mul33x33(real(tensorField_complex(1,1,1,1:3,1:3)),& ! L_2 norm of average stress (freq 0,0,0) in fourier space, pstress_av_L2 = sqrt(maxval (math_eigenvalues33(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress
math_transpose33(real(tensorField_complex(1,1,1,1:3,1:3))))))) ! ignore imaginary part as it is always zero for real only input math_transpose33(pstress_av_lab)))))
err_div_RMS = 0.0_pReal err_div_RMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2) do k = 1_pInt, res(3); do j = 1_pInt, res(2)
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
err_div_RMS = err_div_RMS & err_div_RMS = err_div_RMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
xi(1:3,i,j,k))*two_pi_img)**2.0_pReal)& ! --> sum squared L_2 norm of vector xi(1:3,i,j,k))*two_pi_img)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& +sum(aimag(math_mul33x3_complex(tensorField_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*two_pi_img)**2.0_pReal)) xi(1:3,i,j,k))*two_pi_img)**2.0_pReal))
enddo enddo
err_div_RMS = err_div_RMS & ! Those two layers do not have a conjugate complex counterpart err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart
+ sum(real(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),& + sum(real(math_mul33x3_complex(tensorField_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)& xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),& + sum(aimag(math_mul33x3_complex(tensorField_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)& xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)&
+ sum(real(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),& + sum(real(math_mul33x3_complex(tensorField_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal)& xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),& + sum(aimag(math_mul33x3_complex(tensorField_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal) xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal)
enddo; enddo enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
! if(err_div_RMS/p_hat_avg/sqrt(wgt) * correctionFactor>err_div& if(err_div_RMS/pstress_av_L2*sqrt(wgt) * correctionFactor>err_div&
! .and. iter >2_pInt& .and.iter >2_pInt&
! .and.err_stress > err_stress_tol) iter = itmax .and.err_stress < err_stress_tol) then
err_div = err_div_RMS/p_hat_avg/sqrt(wgt) * correctionFactor ! criterion to stop iterations print*, 'Increasing divergence, stopping iterations'
iter = itmax
endif
err_div = err_div_RMS/pstress_av_L2*sqrt(wgt) * correctionFactor ! criterion to stop iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report ! calculate additional divergence criteria and report
if(debugDivergence) then ! calculate divergence again if(debugDivergence) then ! calculate divergence again
err_div_max = 0.0_pReal err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
temp3_Complex = math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& temp3_Complex = math_mul33x3_complex(tensorField_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*two_pi_img xi(1:3,i,j,k))*two_pi_img
err_div_max = max(err_div_max,sqrt(sum(abs(temp3_Complex)**2.0_pReal))) err_div_max = max(err_div_max,sqrt(sum(abs(temp3_Complex)**2.0_pReal)))
divergence_complex(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared
enddo; enddo; enddo enddo; enddo; enddo
call fftw_execute_dft_c2r(plan_divergence,divergence_complex,divergence_real) call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real)
divergence_real = divergence_real*wgt divergence_real = divergence_real*wgt
err_real_div_RMS = 0.0_pReal err_real_div_RMS = 0.0_pReal
err_real_div_max = 0.0_pReal err_real_div_max = 0.0_pReal
@ -949,8 +943,6 @@ program DAMASK_spectral
err_real_div_RMS = err_real_div_RMS + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of L_2 norm of div(stress) in real space err_real_div_RMS = err_real_div_RMS + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of L_2 norm of div(stress) in real space
err_real_div_max = max(err_real_div_max, sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal))) ! maximum of L two norm of div(stress) in real space err_real_div_max = max(err_real_div_max, sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal))) ! maximum of L two norm of div(stress) in real space
enddo; enddo; enddo enddo; enddo; enddo
p_real_avg = sqrt(maxval (math_eigenvalues33(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space,
math_transpose33(pstress_av_lab)))))
err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space
err_div_max = err_div_max*sqrt(wgt) err_div_max = err_div_max*sqrt(wgt)
@ -960,69 +952,100 @@ program DAMASK_spectral
print '(a,es10.4)', 'error divergence Real max = ',err_real_div_max print '(a,es10.4)', 'error divergence Real max = ',err_real_div_max
print '(a,es10.4)', 'divergence RMS FT/real = ',err_div_RMS/err_real_div_RMS print '(a,es10.4)', 'divergence RMS FT/real = ',err_div_RMS/err_real_div_RMS
print '(a,es10.4)', 'divergence max FT/real = ',err_div_max/err_real_div_max print '(a,es10.4)', 'divergence max FT/real = ',err_div_max/err_real_div_max
print '(a,es10.4)', 'avg stress FT/real = ',p_hat_avg*wgt/p_real_avg
print '(a,es10.4)', 'max deviat. from postProc = ',max_div_error print '(a,es10.4)', 'max deviat. from postProc = ',max_div_error
endif endif
print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', rel. error = ', err_div/err_div_tol print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', rel. error = ', err_div/err_div_tol
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! divergence is calculated from FT(stress), depending on algorithm use field for spectral method ! divergence is calculated from FT(stress), depending on algorithm use field for spectral method
if (.not. simplified_algorithm) tensorField_complex = tau_complex if (.not. simplified_algorithm) tensorField_fourier = tau_fourier
max_diag = tiny(1.0_pReal)
max_offdiag = tiny(1.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium) ! to the actual spectral method calculation (mechanical equilibrium)
if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat
do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red
if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then
do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt do m = 1_pInt ,3_pInt; do p = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) xiDyad(m,p) = xi(m, i,j,k)*xi(p, i,j,k)
enddo; enddo enddo; enddo
temp33_Real = math_inv33(math_mul3333xx33(c0_reference, xiDyad)) ! do l = 1_pInt,3_pInt
! do n = 1_pInt,3_pInt
! temp33_Real(l,n) = sum(c0_reference(l,1:3,n,1:3)*xiDyad)
! enddo; enddo
temp33_Real= math_mul3333xx33(c0_reference,xiDyad)
! max_diag = max(max_diag,maxval( math_mul33x33(temp33_Real,math_inv33(temp33_Real)),&
! reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3])))
! max_offdiag = max(max_offdiag,maxval( math_mul33x33(temp33_Real,math_inv33(temp33_Real)),&
! reshape([ .true.,.false.,.false.,.false.,.true.,.false.,.false.,.false.,.true.],[ 3,3])))
! temp33_Real = math_inv33(temp33_Real)
else else
xiDyad = 0.0_pReal xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
endif endif
do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt
gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*& gamma_hat(1,1,1, l,m,n,p) = 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(xiDyad(m,p) +xiDyad(p,m)) (xiDyad(m,p)+xiDyad(p,m))
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
! if (k==res(3)/2 .or. k==res(3)/2+2 .or.&
! j==res(2)/2 .or. j==res(2)/2+2 .or.&
! i==res(1)/2 .or. i==res(1)/2+2) then
! print*,'gamma_hat', gamma_hat(1,1,1,1:3,1:3,1:3,1:3)
! print*, 's0', s0_reference
! pause
! endif
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
temp33_Complex(m,n) = sum(gamma_hat(1,1,1, m,n, 1:3,1:3) * tensorField_complex(i,j,k,1:3,1:3)) temp33_Complex(m,n) = sum(gamma_hat(1,1,1, m,n, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3))
enddo; enddo enddo; enddo
tensorField_complex(i,j,k,1:3,1:3) = temp33_Complex tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex
enddo; enddo; enddo enddo; enddo; enddo
else ! use precalculated gamma-operator else ! use precalculated gamma-operator
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * tensorField_complex(i,j,k,1:3,1:3)) temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3))
enddo; enddo enddo; enddo
tensorField_complex(i,j,k,1:3,1:3) = temp33_Complex tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex
enddo; enddo; enddo enddo; enddo; enddo
endif endif
tensorField_complex(1,1,1,1:3,1:3) = (defgradAim_lab - defgrad_av_lab)& ! assign average deformation gradient change to zero frequency (real part) tensorField_fourier(1,1,1,1:3,1:3) = (defgrad_av_lab -defgradAim_lab)& ! assign (negative) average deformation gradient change to zero frequency (real part)
* real(Npoints,pReal) * real(Npoints,pReal)
if (.not. simplified_algorithm) tensorField_complex(1,1,1,1:3,1:3) = & ! assign deformation aim to zero frequency (real part) if (.not. simplified_algorithm) tensorField_fourier(1,1,1,1:3,1:3) = & ! assign deformation aim to zero frequency (real part)
defgradAim_lab * real(Npoints,pReal) defgradAim_lab * real(Npoints,pReal)
! print*, 'max off diagonal', max_offdiag
! print*, 'max diagonal', max_diag
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then if (debugFFTW) then
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
scalarField_complex(i,j,k) = tensorField_complex(i,j,k,row,column) scalarField_fourier(i,j,k) = tensorField_fourier(i,j,k,row,column)
enddo; enddo; enddo enddo; enddo; enddo
do i = 0_pInt, res(1)/2_pInt-2_pInt !unpack fft data for conj complex symmetric part. can be directly used in calculation of cstress_field
m = 1_pInt
do k = 1_pInt, res(3)
n = 1_pInt
do j = 1_pInt, res(2)
scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m))
if(n == 1_pInt) n = res(2) + 1_pInt
n = n-1_pInt
enddo
if(m == 1_pInt) m = res(3) + 1_pInt
m = m -1_pInt
enddo; enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing the inverse FT ! doing the inverse FT
call fftw_execute_dft_c2r(plan_correction,tensorField_complex,tensorField_real) ! back transform of fluct deformation gradient call fftw_execute_dft_c2r(plan_correction,tensorField_fourier,tensorField_real) ! back transform of fluct deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results ! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then if (debugFFTW) then
print '(a,i1,1x,i1)', 'checking iFT results of compontent ', row, column print '(a,i1,1x,i1)', 'checking iFT results of compontent ', row, column
call fftw_execute_dft_c2r(plan_scalarField_back,scalarField_complex,scalarField_real) call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real)
print '(a,es10.4)', 'max iFT relative error ',& print '(a,es10.4)', 'max iFT relative error ',&
maxval((scalarField_real(1:res(1),1:res(2),1:res(3))-& maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-&
tensorField_real(1:res(1),1:res(2),1:res(3),row,column))/& tensorField_real(1:res(1),1:res(2),1:res(3),row,column))/&
scalarField_real(1:res(1),1:res(2),1:res(3))) real(scalarField_real(1:res(1),1:res(2),1:res(3))))
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1049,7 +1072,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient ! updated deformation gradient
defgrad = defgrad + tensorField_real(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization defgrad = defgrad - tensorField_real(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient in case of fluctuation field algorithm ! updated deformation gradient in case of fluctuation field algorithm
@ -1080,16 +1103,37 @@ program DAMASK_spectral
print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged' print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
else else
convergedCounter = convergedCounter + 1_pInt
print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' converged' print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' converged'
endif endif
if (mod(totalIncsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency if (mod(totalIncsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency
print '(a)', '' print '(a)', ''
print '(a)', '... writing results to file ......................................' print '(a)', '... writing results to file ......................................'
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
endif endif
if( bc(loadcase)%restartFrequency > 0_pInt .and. &
mod(inc - 1_pInt,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartWrite = .true.
print '(A)', 'writing converged results for restart'
if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! writing deformation gradient field to file
write (777,rec=1) defgrad
close (777)
endif
restartInc=totalIncsCounter
endif
if (update_gamma) then if (update_gamma) then
print*, 'update c0_reference ' print*, 'update c0_reference '
c0_reference = c_current*wgt c0_reference = c_current*wgt
! s0_reference = math_Plain99to3333(s0_99)
!c0_99 = math_Plain3333to99(c0_reference)
! call math_invert(9_pInt, s0_99, c0_99, i, errmatinv) ! invert reduced stiffness
! if(errmatinv) call IO_error(error_ID=800_pInt)
! print*, (c0_reference - math_Plain99to3333(c0_99))/c0_reference
! pause
endif endif
endif ! end calculation/forwarding endif ! end calculation/forwarding
@ -1100,7 +1144,7 @@ program DAMASK_spectral
print '(a)', '' print '(a)', ''
print '(a)', '##################################################################' print '(a)', '##################################################################'
print '(i6.6,a,i6.6,a)', notConvergedCounter, ' out of ', & print '(i6.6,a,i6.6,a)', notConvergedCounter, ' out of ', &
totalIncsCounter - restartReadInc, ' increments did not converge!' notConvergedCounter + convergedCounter, ' increments did not converge!'
close(538) close(538)
call fftw_destroy_plan(plan_stress); call fftw_destroy_plan(plan_correction) call fftw_destroy_plan(plan_stress); call fftw_destroy_plan(plan_correction)
if (debugDivergence) call fftw_destroy_plan(plan_divergence) if (debugDivergence) call fftw_destroy_plan(plan_divergence)

42
code/FEsolving.f90
View File

@ -21,11 +21,10 @@
!############################################################## !##############################################################
MODULE FEsolving MODULE FEsolving
!############################################################## !##############################################################
use prec, only: pInt,pReal use prec, only: pInt,pReal
implicit none implicit none
integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartReadInc = 0_pInt integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartInc = 1_pInt
real(pReal) :: theTime = 0.0_pReal, theDelta = 0.0_pReal real(pReal) :: theTime = 0.0_pReal, theDelta = 0.0_pReal
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false. logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false.
logical :: symmetricSolver = .false. logical :: symmetricSolver = .false.
@ -46,6 +45,7 @@
!*********************************************************** !***********************************************************
subroutine FE_init() subroutine FE_init()
use, intrinsic :: iso_fortran_env
use prec, only: pInt use prec, only: pInt
use debug, only: debug_verbosity use debug, only: debug_verbosity
use DAMASK_interface use DAMASK_interface
@ -55,7 +55,7 @@
integer(pInt), parameter :: fileunit = 222 integer(pInt), parameter :: fileunit = 222
integer(pInt), parameter :: maxNchunks = 6 integer(pInt), parameter :: maxNchunks = 6
integer(pInt):: i, start = 0_pInt, length=0_pInt integer(pInt):: i, start = 0_pInt, length=0_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
character(len=64) tag character(len=64) tag
character(len=1024) line, commandLine character(len=1024) line, commandLine
@ -74,27 +74,29 @@
if(start /= 0_pInt) then ! found something if(start /= 0_pInt) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartReadInc ! read argument read(commandLine(start:start+length),'(I12)') restartInc ! read argument
restartReadInc = restartReadInc - 1_pInt ! command line argument is inc to compute restartRead = restartInc > 0_pInt
restartRead = max(0_pInt,restartReadInc) > 0_pInt if(restartInc <= 0_pInt) then
if(restartReadInc < 0_pInt) call IO_warning(warning_ID=34_pInt) call IO_warning(warning_ID=34_pInt)
restartInc = 1_pInt
endif
endif endif
else else
rewind(fileunit) rewind(fileunit)
do do
read (fileunit,'(a1024)',END=100) line read (fileunit,'(a1024)',END=100) line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag) select case(tag)
case ('solver') case ('solver')
read (fileunit,'(a1024)',END=100) line ! next line read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt) symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt)
case ('restart') case ('restart')
read (fileunit,'(a1024)',END=100) line ! next line read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
restartWrite = iand(IO_intValue(line,positions,1),1_pInt) > 0_pInt restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1),2_pInt) > 0_pInt restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart') case ('*restart')
do i=2,positions(1) do i=2,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite
@ -109,7 +111,7 @@
enddo enddo
endif endif
else else
call IO_error(101, ext_msg=FEmodelGeometry) ! cannot open input file call IO_error(101_pInt, ext_msg=FEmodelGeometry) ! cannot open input file
endif endif
100 close(fileunit) 100 close(fileunit)
@ -119,27 +121,27 @@
do do
read (fileunit,'(a1024)',END=200) line read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,positions,1)) == 'restart' .and. & if ( IO_lc(IO_stringValue(line,positions,1_pInt)) == 'restart' .and. &
IO_lc(IO_stringValue(line,positions,2)) == 'file' .and. & IO_lc(IO_stringValue(line,positions,2_pInt)) == 'file' .and. &
IO_lc(IO_stringValue(line,positions,3)) == 'job' .and. & IO_lc(IO_stringValue(line,positions,3_pInt)) == 'job' .and. &
IO_lc(IO_stringValue(line,positions,4)) == 'id' ) & IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) &
FEmodelGeometry = IO_StringValue(line,positions,6) FEmodelGeometry = IO_StringValue(line,positions,6_pInt)
enddo enddo
elseif (FEsolver == 'Abaqus' .and. IO_open_inputFile(fileunit,FEmodelGeometry)) then elseif (FEsolver == 'Abaqus' .and. IO_open_inputFile(fileunit,FEmodelGeometry)) then
rewind(fileunit) rewind(fileunit)
do do
read (fileunit,'(a1024)',END=200) line read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,positions,1))=='*heading') then if ( IO_lc(IO_stringValue(line,positions,1_pInt))=='*heading') then
read (fileunit,'(a1024)',END=200) line read (fileunit,'(a1024)',END=200) line
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
FEmodelGeometry = IO_StringValue(line,positions,1) FEmodelGeometry = IO_StringValue(line,positions,1_pInt)
endif endif
enddo enddo
elseif (FEsolver == 'Spectral') then elseif (FEsolver == 'Spectral') then
!do nothing !do nothing
else else
call IO_error(106) ! cannot open file for old job info call IO_error(106_pInt) ! cannot open file for old job info
endif endif
endif endif

133
code/makefile
View File

@ -33,8 +33,9 @@
# SUFFIX = arbitrary suffix # SUFFIX = arbitrary suffix
# STANDARD_CHECK = checking for Fortran 2008, compiler dependend # STANDARD_CHECK = checking for Fortran 2008, compiler dependend
######################################################################################## ########################################################################################
# Here are some useful debugging switches. Switch on by uncommenting the #SUFFIX line at the end of this section: # Here are some useful debugging switches for ifort. Switch on by uncommenting the #SUFFIX line at the end of this section:
# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ # information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/
# check if an array index is too small (<1) or too large! # check if an array index is too small (<1) or too large!
DEBUG1 =-check bounds -g DEBUG1 =-check bounds -g
@ -54,9 +55,12 @@ DEBUG5 =-warn all
#set precision (check for missing _pInt and _pReal) #set precision (check for missing _pInt and _pReal)
DEBUG6= -real-size 32 -integer-size 16 DEBUG6= -real-size 32 -integer-size 16
#or one of those 16/32/64/128 (= 2,4,8,16 bytes) #or one of those 16/32/64/128 (= 2,4,8,16 bytes)
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3) $(DEBUG4) $(DEBUG5) $(DEBUG6)
# Here are some useful debugging switches for gfortran
# fcheck-bounds: eqv to DEBUG1 of ifort
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3)
######################################################################################## ########################################################################################
#auto values will be set by setup_code.py #auto values will be set by setup_code.py
@ -70,6 +74,12 @@ COMPILERNAME ?= $(F90)
OPENMP ?= ON OPENMP ?= ON
OPTIMIZATION ?= DEFENSIVE OPTIMIZATION ?= DEFENSIVE
ifeq "$(F90)" "ifort"
ARCHIVE_COMMAND :=xiar
else
ARCHIVE_COMMAND :=ar
endif
ifeq "$(OPTIMIZATION)" "OFF" ifeq "$(OPTIMIZATION)" "OFF"
OPTI := OFF OPTI := OFF
MAXOPTI := OFF MAXOPTI := OFF
@ -93,40 +103,40 @@ MAXOPTI := DEFENSIVE
endif endif
ifeq "$(PORTABLE)" "FALSE" ifeq "$(PORTABLE)" "FALSE"
PORTABLE_SWITCH = -msse3 PORTABLE_SWITCH =-msse3
endif endif
# settings for multicore support # settings for multicore support
ifeq "$(OPENMP)" "ON" ifeq "$(OPENMP)" "ON"
OPENMP_FLAG_ifort = -openmp -openmp-report0 -parallel OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel
OPENMP_FLAG_gfortran = -fopenmp OPENMP_FLAG_gfortran =-fopenmp
ACML_ARCH =_mp ACML_ARCH =_mp
LIBRARIES += -lfftw3_threads -lpthread LIBRARIES +=-lfftw3_threads -lpthread
endif endif
LIBRARIES += -lfftw3 LIBRARIES +=-lfftw3
LIB_DIRS += -L$(FFTWROOT)/lib LIB_DIRS +=-L$(FFTWROOT)/lib
ifdef IKMLROOT ifdef IKMLROOT
LIBRARIES += -mkl LIBRARIES +=-mkl
else else
ifdef ACMLROOT ifdef ACMLROOT
LIB_DIRS += -L$(ACMLROOT)/$(F90)64$(ACML_ARCH)/lib LIB_DIRS +=-L$(ACMLROOT)/$(F90)64$(ACML_ARCH)/lib
LIBRARIES += -lacml$(ACML_ARCH) LIBRARIES +=-lacml$(ACML_ARCH)
else else
ifdef LAPACKROOT ifdef LAPACKROOT
LIB_DIRS += -L$(LAPACKROOT)/lib64 -L$(LAPACKROOT)/lib LIB_DIRS +=-L$(LAPACKROOT)/lib64 -L$(LAPACKROOT)/lib
LIBRARIES += -llapack LIBRARIES +=-llapack
endif endif
endif endif
endif endif
ifdef STANDARD_CHECK ifdef STANDARD_CHECK
STANDARD_CHECK_ifort = $(STANDARD_CHECK) STANDARD_CHECK_ifort =$(STANDARD_CHECK)
STANDARD_CHECK_gfortran = $(STANDARD_CHECK) STANDARD_CHECK_gfortran =$(STANDARD_CHECK)
endif endif
STANDARD_CHECK_ifort ?= -stand f08 STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics
STANDARD_CHECK_gfortran ?=-std=f2008 STANDARD_CHECK_gfortran ?=-std=f2008
@ -138,16 +148,89 @@ OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ip -static -fp-model
OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize
COMPILE_OPTIONS_ifort := -fpp -diag-disable 8291,8290,5268 COMPILE_OPTIONS_ifort :=-fpp\
#warning ID 9291,8290: -diag-enable sc3\
#warning ID 5268: Extension to standard: The text exceeds right hand column allowed on the line (we have only comments there) -diag-disable 8291,8290,5268\
COMPILE_OPTIONS_gfortran := -xf95-cpp-input -ffree-line-length-132 -fno-range-check -warn declarations\
-warn general\
-warn usage
#alignments: Determines whether warnings occur for data that is not naturally aligned.
#declarations: Determines whether warnings occur for any undeclared names.
#errors: Determines whether warnings are changed to errors.
#general: Determines whether warning messages and informational messages are issued by the compiler.
#ignore_loc: Determines whether warnings occur when %LOC is stripped from an actual argument.
#interfaces: Determines whether the compiler checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks.
#stderrors: Determines whether warnings about Fortran standard violations are changed to errors.
#truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files.
#uncalled: Determines whether warnings occur when a statement function is never called
#unused: Determines whether warnings occur for declared variables that are never used.
#usage: Determines whether warnings occur for questionable programming practices.
#-fpp: preprocessor
#-diag-disable: disables warnings, where
# warning ID 9291:
# warning ID 8290:
# warning ID 5268: The text exceeds right hand column allowed on the line (we have only comments there)
COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-ffree-line-length-132\
-fno-range-check\
-fimplicit-none\
-pedantic\
-Warray-bounds\
-Wunused-parameter\
-Wampersand\
-Wno-tabs\
-Wcharacter-truncation\
-Wintrinsic-shadow\
-Waliasing\
-Wconversion\
-Wsurprising\
-Wunused-value\
-Wunderflow
#-xf95-cpp-input: preprocessor
#-ffree-line-length-132: restrict line length to the standard 132 characters
#-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN
#-fimplicit-none: assume "implicit-none" even if not present in source
#-pedantic: more strict on standard, enables some of the warnings below
#-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime
#-Wunused-parameter: find usused variables with "parameter" attribute
#-Wampersand: checks if a character expression is continued proberly by an ampersand at the end of the line and at the beginning of the new line
#-Wno-tabs: do not allow tabs in source
#-Wcharacter-truncation: warn if character expressions (strings) are truncated
#-Wintrinsic-shadow: warn if a user-defined procedure or module procedure has the same name as an intrinsic
#-Waliasing: warn about possible aliasing of dummy arguments. Specifically, it warns if the same actual argument is associated with a dummy argument with "INTENT(IN)" and a dummy argument with "INTENT(OUT)" in a call with an explicit interface.
#-Wconversion: warn about implicit conversions between different type
#-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made.
#-Wunused-value:
#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation
#MORE OPTIONS
# only for gfortran 4.6:
#-Wsuggest-attribute=const
#-Wsuggest-attribute=noreturn
#-Wsuggest-attribute=pure
# too many warnings because we have comments beyond character 132:
#-Wline-truncation
# warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers:
#-Wintrinsic-std
# warnings because we have many temporary arrays (performance issue?):
#-Warray-temporaries
# -Wimplicit-interface
# -pedantic-errors
# -fmodule-private
COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c
COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c
DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a
$(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a \ $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a \
constitutive.a advanced.a basics.a $(LIB_DIRS) $(LIBRARIES) constitutive.a advanced.a basics.a $(LIB_DIRS) $(LIBRARIES)
@ -158,7 +241,7 @@ DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o
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 $(ARCHIVE_COMMAND) rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o
CPFEM.o: CPFEM.f90 homogenization.o CPFEM.o: CPFEM.f90 homogenization.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) CPFEM.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) CPFEM.f90 $(SUFFIX)
@ -178,7 +261,7 @@ crystallite.o: crystallite.f90 constitutive.a
constitutive.a: constitutive.o constitutive.a: constitutive.o
ar rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a $(ARCHIVE_COMMAND) rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a
constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive.f90 $(SUFFIX)
@ -201,7 +284,7 @@ constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a
advanced.a: lattice.o advanced.a: lattice.o
ar rc advanced.a FEsolving.o mesh.o material.o lattice.o $(ARCHIVE_COMMAND) rc advanced.a FEsolving.o mesh.o material.o lattice.o
lattice.o: lattice.f90 material.o lattice.o: lattice.f90 material.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) lattice.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) lattice.f90 $(SUFFIX)
@ -218,7 +301,7 @@ FEsolving.o: FEsolving.f90 basics.a
basics.a: math.o basics.a: math.o
ar rc basics.a math.o debug.o numerics.o IO.o DAMASK_spectral_interface.o prec.o $(ARCHIVE_COMMAND) rc basics.a math.o debug.o numerics.o IO.o DAMASK_spectral_interface.o prec.o
math.o: math.f90 debug.o math.o: math.f90 debug.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) math.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) math.f90 $(SUFFIX)