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

removed phase contrast/preconditioning 

added information on itmin in output, impoved output
set exit code to 0 on successful termination (seems to be unix standard)
exit codes:
0: successful termination
1: error (using IO_error)
2: require regrid 
updated the AL solver, still VERY experimental
This commit is contained in:
Martin Diehl 2012-06-05 16:34:20 +00:00
parent 4b6800b89a
commit fa7f9866df
3 changed files with 166 additions and 190 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
code/DAMASK_spectral.f90
View File

@ -195,7 +195,6 @@ program DAMASK_spectral
real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc
real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature real(pReal), dimension(:,:,:), allocatable :: temperature
real(pReal), dimension(:,:,:), allocatable :: phase_cont ! phase contrast field: C(x)/C_ref
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW ! variables storing information for spectral method and FFTW
@ -449,7 +448,6 @@ program DAMASK_spectral
allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally
allocate (phase_cont ( res(1), res(2),res(3)), source = 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, P_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, P_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, deltaF_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, deltaF_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
@ -626,29 +624,29 @@ C_ref = C * wgt
!################################################################################################## !##################################################################################################
do inc = 1_pInt, bc(loadcase)%incs do inc = 1_pInt, bc(loadcase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1_pInt
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeinc_old = timeinc timeinc_old = timeinc
if (bc(loadcase)%logscale == 0_pInt) then ! linear scale 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 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 else
if (loadcase == 1_pInt) then ! 1st loadcase of 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 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 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 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)) 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
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 endif
time = time + timeinc
if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding)
if (bc(loadcase)%velGradApplied) then ! calculate deltaF_aim from given L and current F if (bc(loadcase)%velGradApplied) then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * mask_defgrad * math_mul33x33(bc(loadcase)%deformation, F_aim) deltaF_aim = timeinc * mask_defgrad * math_mul33x33(bc(loadcase)%deformation, F_aim)
else ! deltaF_aim = fDot *timeinc where applicable else ! deltaF_aim = fDot *timeinc where applicable
@ -729,8 +727,8 @@ C_ref = C * wgt
! report begin of new iteration ! report begin of new iteration
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
write(6,'(5(a,i6.6))') 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,& write(6,'(6(a,i6.6))') 'Loadcase ',loadcase,' Inc. ',inc,'/',bc(loadcase)%incs,&
' @ Iteration ',iter,'/',itmax ' @ Iter. ',itmin,' < ',iter,' < ',itmax
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',&
math_transpose33(F_aim) math_transpose33(F_aim)
write(6,'(a)') '' write(6,'(a)') ''
@ -747,8 +745,6 @@ C_ref = C * wgt
P_real(i,j,k,1:3,1:3),dPdF) P_real(i,j,k,1:3,1:3),dPdF)
enddo; enddo; enddo enddo; enddo; enddo
P_real = 0.0_pReal ! needed because of the padding for FFTW P_real = 0.0_pReal ! needed because of the padding for FFTW
C = 0.0_pReal C = 0.0_pReal
ielem = 0_pInt ielem = 0_pInt
@ -760,11 +756,10 @@ C_ref = C * wgt
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, & temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, &
P_real(i,j,k,1:3,1:3),dPdF) P_real(i,j,k,1:3,1:3),dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
phase_cont(i,j,k) = sqrt(sum(dPdF*dPdF)/sum(C_ref*C_ref))
C = C + dPdF C = C + dPdF
enddo; enddo; enddo enddo; enddo; enddo
call debug_info() call debug_info()
write(6,'(a,es11.4)') 'Max phase contrast = ',maxval(phase_cont)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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
@ -814,8 +809,6 @@ C_ref = C * wgt
P_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& P_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) = 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
@ -864,8 +857,9 @@ C_ref = C * wgt
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/pstress_av_L2 > err_div & if (err_div_RMS/pstress_av_L2 > err_div &
.and. err_stress < err_stress_tol) then .and. err_stress < err_stress_tol &
.and. iter >= itmin ) then
write(6,'(a)') 'Increasing divergence, stopping iterations' write(6,'(a)') 'Increasing divergence, stopping iterations'
iter = itmax iter = itmax
endif endif
@ -1000,8 +994,8 @@ C_ref = C * wgt
! updated deformation gradient ! updated deformation gradient
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)
F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - deltaF_real(i,j,k,1:3,1:3)*wgt/phase_cont(i,j,k) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - deltaF_real(i,j,k,1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
enddo; enddo; enddo ! preconditioning: F(x)^(n+1) = F(x)^(n) + correction/phase_contrast(x) enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1059,8 +1053,11 @@ C_ref = C * wgt
enddo ! end looping over loadcases enddo ! end looping over loadcases
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', & write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' increments did not converge!' notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
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)
@ -1068,5 +1065,5 @@ C_ref = C * wgt
call fftw_destroy_plan(plan_scalarField_forth) call fftw_destroy_plan(plan_scalarField_forth)
call fftw_destroy_plan(plan_scalarField_back) call fftw_destroy_plan(plan_scalarField_back)
endif endif
call quit(1_pInt) call quit(0_pInt)
end program DAMASK_spectral end program DAMASK_spectral

279
code/DAMASK_spectral_AL.f90
View File

@ -89,7 +89,6 @@ program DAMASK_spectral_AL
restartInc restartInc
use numerics, only: & use numerics, only: &
err_div_tol, &
err_stress_tolrel, & err_stress_tolrel, &
err_stress_tolabs, & err_stress_tolabs, &
rotation_tol, & rotation_tol, &
@ -158,11 +157,11 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: P_star_av = 0.0_pReal, & real(pReal), dimension(3,3) :: P_av = 0.0_pReal, &
F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, & F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, &
mask_stress, mask_defgrad, deltaF, F_star_av, & mask_stress, mask_defgrad, deltaF, F_star_av, &
F_aim_lab ! quantities rotated to other coordinate system F_aim_lab ! quantities rotated to other coordinate system
real(pReal), dimension(3,3,3,3) :: C_inc0, C=0.0_pReal, S_lastInc, C_lastInc ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: C_inc0, C=0.0_pReal, S_inc0, S_lastInc, C_lastInc ! stiffness and compliance
real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde ! 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 ! compliance and stiffness in matrix notation
@ -172,16 +171,16 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! pointwise data ! pointwise data
type(C_PTR) :: tensorField ! fields in real an fourier space type(C_PTR) :: tensorField ! fields in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: lambda_real, F_real ! fields in real space (pointer) real(pReal), dimension(:,:,:,:,:), pointer :: tau_real, F_real ! fields in real space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: lambda_fourier, F_fourier ! fields in fourier space (pointer) complex(pReal), dimension(:,:,:,:,:), pointer :: tau_fourier, F_fourier ! fields in fourier space (pointer)
real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, F_star, lambda, P, F_star_lastIter real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, lambda, P, F_star
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: dPdF real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: dPdF
real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature real(pReal), dimension(:,:,:), allocatable :: temperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW ! variables storing information for spectral method and FFTW
type(C_PTR) :: plan_correction, plan_lambda ! plans for fftw type(C_PTR) :: plan_correction, plan_tau ! plans for fftw
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method
real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator
@ -191,7 +190,7 @@ program DAMASK_spectral_AL
! loop variables, convergence etc. ! loop variables, convergence etc.
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) :: 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_stress, err_stress_tol, err_f, err_p, err_crit, & real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit, &
err_f_point, err_p_point, pstress_av_L2, err_div_rms, err_div err_f_point, err_p_point, err_nl, err_nl_point
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,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: temp33_Real real(pReal), dimension(3,3) :: temp33_Real
@ -325,6 +324,8 @@ program DAMASK_spectral_AL
res = mesh_spectral_getResolution() res = mesh_spectral_getResolution()
geomdim = mesh_spectral_getDimension() geomdim = mesh_spectral_getDimension()
homog = mesh_spectral_getHomogenization() homog = mesh_spectral_getHomogenization()
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
Npoints = res(1)*res(2)*res(3) Npoints = res(1)*res(2)*res(3)
wgt = 1.0_pReal/real(Npoints, pReal) wgt = 1.0_pReal/real(Npoints, pReal)
@ -407,16 +408,15 @@ program DAMASK_spectral_AL
allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (dPdF ( res(1), res(2),res(3),3,3,3,3), source = 0.0_pReal) allocate (dPdF ( res(1), res(2),res(3),3,3,3,3), source = 0.0_pReal)
allocate (F_star ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_star ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (F_star_lastIter ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally
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, lambda_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation call c_f_pointer(tensorField, tau_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation
call c_f_pointer(tensorField, F_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation call c_f_pointer(tensorField, F_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation
call c_f_pointer(tensorField, lambda_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation call c_f_pointer(tensorField, tau_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation
call c_f_pointer(tensorField, F_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation call c_f_pointer(tensorField, F_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -431,10 +431,10 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! creating plans ! creating plans
plan_lambda = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order plan_tau = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order
lambda_real,[ res(3),res(2) ,res(1)+2_pInt],& ! input data , physical length in each dimension in reversed order tau_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 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
lambda_fourier,[ res(3),res(2) ,res1_red],& tau_fourier,[ res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,fftw_planner_flag) 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,& plan_correction = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],9,&
@ -449,7 +449,8 @@ program DAMASK_spectral_AL
ielem = 0_pInt ielem = 0_pInt
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 ielem = ielem + 1_pInt
F_real(i,j,k,1:3,1:3) = math_I3; F_lastInc(i,j,k,1:3,1:3) = math_I3 F_star(i,j,k,1:3,1:3) = math_I3
F_lastInc(i,j,k,1:3,1:3) = math_I3
coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal) coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal)
call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3)) 0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3))
@ -462,7 +463,7 @@ program DAMASK_spectral_AL
0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3)) 0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3))
C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3) C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
C_inc0 = C * wgt ! linear reference material stiffness C_inc0 = C * wgt *3.0_pReal ! 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
@ -497,24 +498,6 @@ program DAMASK_spectral_AL
gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
endif endif
!--------------------------------------------------------------------------------------------------
! possible restore deformation gradient from saved state
if (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F_star))
read (777,rec=1) F_star
close (777)
F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_star
F_lastInc = F_star
F_aim = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
F_aim = F_aim + F_real(i,j,k,1:3,1:3) ! calculating old average deformation
enddo; enddo; enddo
F_aim = F_aim * wgt
F_aim_lastInc = F_aim
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
@ -601,15 +584,14 @@ program DAMASK_spectral_AL
+ guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & + guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old &
+ deltaF + deltaF
F_aim_lastInc = temp33_Real F_aim_lastInc = temp33_Real
F_star_av = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update local deformation gradient ! update local deformation gradient
deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation) 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 = F_real(i,j,k,1:3,1:3) temp33_Real = F_star(i,j,k,1:3,1:3)
F_real(i,j,k,1:3,1:3) = F_real(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (F_real(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))& ! guessing... + guessmode * (F_star(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))& ! guessing...
*timeinc/timeinc_old & *timeinc/timeinc_old &
+ (1.0_pReal-guessmode) * deltaF ! if not guessing, use prescribed average deformation where applicable + (1.0_pReal-guessmode) * deltaF ! if not guessing, use prescribed average deformation where applicable
F_lastInc(i,j,k,1:3,1:3) = temp33_Real F_lastInc(i,j,k,1:3,1:3) = temp33_Real
@ -617,19 +599,23 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Initialize / Update lambda to useful value ! Initialize / Update lambda to useful value
P_star_av = P_star_av + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc) P_av = P_av + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc)
lambda_av = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
lambda(i,j,k,1:3,1:3) = P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3), &
F_real(i,j,k,1:3,1:3)-F_lastInc(i,j,k,1:3,1:3))
lambda_av = lambda_av + lambda(i,j,k,1:3,1:3)
enddo; enddo; enddo
lambda_av=lambda_av*wgt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc)
enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!Initialize pointwise data for AL scheme: ToDo: good choice? !Initialize pointwise data for AL scheme: ToDo: good choice?
F_star(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3)
F_star_av = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
F_star_av = F_star_av + F_star(i,j,k,1:3,1:3)
enddo; enddo; enddo
F_star_av = F_star_av *wgt
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',&
math_transpose33(F_star_av)
!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate reduced compliance ! calculate reduced compliance
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
@ -660,7 +646,7 @@ program DAMASK_spectral_AL
endif; enddo; endif; enddo endif; enddo; endif; enddo
S_lastInc = (math_Plain99to3333(s_prev99)) S_lastInc = (math_Plain99to3333(s_prev99))
endif endif
S_inc0 = math_invSym3333(C*wgt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
write(6,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
@ -690,9 +676,9 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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
err_stress = maxval(abs(mask_stress * (P_star_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable)
F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- P_star_av) F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- P_av)
err_stress_tol = maxval(abs(P_star_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
else else
err_stress_tol = + huge(1.0_pReal) err_stress_tol = + huge(1.0_pReal)
endif endif
@ -703,42 +689,19 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing Fourier transform ! doing Fourier transform
write(6,'(a)') '... spectral method ...............................................' write(6,'(a)') '... spectral method ...............................................'
lambda_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = lambda(1:res(1),1:res(2),1:res(3),1:3,1:3) tau_real = 0.0_pReal
call fftw_execute_dft_r2c(plan_lambda,lambda_real,lambda_fourier) do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res(1)
lambda_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& tau_real(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) - math_mul3333xx33(C_inc0,F_star(i,j,k,1:3,1:3))
enddo; enddo; enddo
call fftw_execute_dft_r2c(plan_tau,tau_real,tau_fourier)
tau_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
lambda_fourier(1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& tau_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) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)>1_pInt) & if(res(3)>1_pInt) &
lambda_fourier(1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& tau_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) = cmplx(0.0_pReal,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
pstress_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(lambda_av,& ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
math_transpose33(lambda_av)))))
err_div_RMS = 0.0_pReal
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.
err_div_RMS = err_div_RMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(lambda_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))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(lambda_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo
err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart
+ sum( real(math_mul33x3_complex(lambda_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(lambda_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum( real(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt
err_div = err_div_RMS/pstress_av_L2
!--------------------------------------------------------------------------------------------------
! using gamma operator to update F ! using gamma operator to update F
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
@ -752,7 +715,7 @@ program DAMASK_spectral_AL
gamma_hat(1,1,1, l,u,v,w) = temp33_Real(l,v)*xiDyad(u,w) gamma_hat(1,1,1, l,u,v,w) = temp33_Real(l,v)*xiDyad(u,w)
forall(l = 1_pInt:3_pInt, u = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, u = 1_pInt:3_pInt) &
temp33_Complex(l,u) = sum(gamma_hat(1,1,1, l,u, 1:3,1:3) *& temp33_Complex(l,u) = sum(gamma_hat(1,1,1, l,u, 1:3,1:3) *&
lambda_fourier(i,j,k,1:3,1:3)) tau_fourier(i,j,k,1:3,1:3))
F_fourier(i,j,k,1:3,1:3) = - temp33_Complex F_fourier(i,j,k,1:3,1:3) = - temp33_Complex
endif endif
enddo; enddo; enddo enddo; enddo; enddo
@ -760,24 +723,39 @@ program DAMASK_spectral_AL
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
forall( u = 1_pInt:3_pInt, v = 1_pInt:3_pInt) & forall( u = 1_pInt:3_pInt, v = 1_pInt:3_pInt) &
temp33_Complex(u,v) = sum(gamma_hat(i,j,k, u,v, 1:3,1:3) *& temp33_Complex(u,v) = sum(gamma_hat(i,j,k, u,v, 1:3,1:3) *&
lambda_fourier(i,j,k,1:3,1:3)) tau_fourier(i,j,k,1:3,1:3))
F_fourier(i,j,k, 1:3,1:3) = - temp33_Complex F_fourier(i,j,k, 1:3,1:3) = - temp33_Complex
enddo; enddo; enddo enddo; enddo; enddo
endif endif
F_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab - F_star_av)*real(Npoints,pReal),0.0_pReal,pReal) F_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab)*real(Npoints,pReal),0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing inverse Fourier transform ! doing inverse Fourier transform
call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real) ! back transform of fluct deformation gradient call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real) ! back transform of fluct deformation gradient
F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) * wgt + & F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) * wgt
F_star(1:res(1),1:res(2),1:res(3),1:3,1:3)
F_star_av = 0.0_pReal
temp33_real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
F_star_av = F_star_av + F_star(i,j,k,1:3,1:3)
temp33_real = temp33_real + F_real(i,j,k,1:3,1:3)
enddo; enddo; enddo
F_star_av = F_star_av *wgt
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',&
math_transpose33(F_star_av)
temp33_real = temp33_real *wgt
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F =',&
math_transpose33(temp33_real)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................'
if(callCPFEM) then err_nl_point = huge(1.0_pReal)
write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................' u = 0_pInt
F_star_lastIter = F_star do while (err_nl_point>1.0e6_pReal .or. u<2_pInt)
u = u + 1_pInt
ielem = 0_pInt ielem = 0_pInt
err_nl_point = 0.0_pReal
err_nl = 0.0_pReal
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 ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle call CPFEM_general(3_pInt,& ! collect cycle
@ -793,83 +771,93 @@ program DAMASK_spectral_AL
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3)) sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3))
CPFEM_mode = 2_pInt ! winding forward CPFEM_mode = 2_pInt ! winding forward
temp33_Real = lambda(i,j,k,1:3,1:3) - P(i,j,k,1:3,1:3) & temp33_Real = P(i,j,k,1:3,1:3) - lambda(i,j,k,1:3,1:3) &
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3)) + math_mul3333xx33(C_inc0,F_star(i,j,k,1:3,1:3)-F_real(i,j,k,1:3,1:3))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + math_mul3333xx33(math_invSym3333(& err_nl_point = max(err_nl_point, maxval(abs(temp33_real)))
C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real) err_nl = max(err_nl, sqrt(math_mul33xx33(temp33_real,temp33_real)))
! temp33_Real = lambda(i,j,k,1:3,1:3) - P(i,j,k,1:3,1:3) &
! + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
! F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + math_mul3333xx33(math_invSym3333(&
! C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real)
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) - 0.1_pReal * math_mul3333xx33(S_inc0, temp33_Real)
enddo; enddo; enddo enddo; enddo; enddo
else
guesses = guesses +1_pInt if(u>1) print*, 'comp wise error of nl eq last it:', err_nl_point/1.0e6_pReal
write(6,'(a)')' ... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax if(u>1) print*, 'norm error: last it ', err_nl/1.0e6_pReal
enddo
write(6,'(a)') '... update λ..........................'
err_f = 0.0_pReal
err_f_point = 0.0_pReal
err_p = 0.0_pReal
err_p_point = 0.0_pReal
F_star_av = 0.0_pReal
lambda_av = 0.0_pReal
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 = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) &
F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))& - F_star(i,j,k,1:3,1:3))
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3)) F_star_av = F_star_av + F_star(i,j,k,1:3,1:3)
lambda_av = lambda_av + lambda(i,j,k,1:3,1:3)
temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3)
err_f_point = max(err_f_point, maxval(abs(temp33_real)))
err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real)))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + math_mul3333xx33(math_invSym3333(&
C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real)
enddo; enddo; enddo enddo; enddo; enddo
endif
write(6,'(a)') '... update λ..........................' ielem = 0_pInt
P_av = 0.0_pReal
err_f = 0.0_pReal C = 0.0_pReal
err_f_point = 0.0_pReal
err_p = 0.0_pReal
err_p_point = 0.0_pReal
F_star_av = 0.0_pReal
P_star_av = 0.0_pReal
lambda_av = 0.0_pReal
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)
lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) & ielem = ielem + 1_pInt
- F_star(i,j,k,1:3,1:3)) call CPFEM_general(3_pInt,& ! collect cycle
F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
lambda_av = lambda_av + lambda(i,j,k,1:3,1:3) F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
P_star_av = P_star_av + P(i,j,k,1:3,1:3) sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3))
enddo; enddo; enddo
temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3) ielem = 0_pInt
err_f_point = max(err_f_point, maxval(abs(temp33_real))) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real))) ielem = ielem + 1_pInt
call CPFEM_general(2_pInt,&
temp33_real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3))) F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3))
P_av = P_av + P(i,j,k,1:3,1:3)
temp33_real = lambda(i,j,k,1:3,1:3) - P(i,j,k,1:3,1:3)
err_p_point = max(err_p_point, maxval(abs(temp33_real))) err_p_point = max(err_p_point, maxval(abs(temp33_real)))
err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real))) err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real)))
C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
F_star_av = F_star_av *wgt F_star_av = F_star_av *wgt
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',&
math_transpose33(F_star_av) math_transpose33(F_star_av)
P_star_av = P_star_av *wgt P_av = P_av *wgt
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) / GPa =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) / GPa =',&
math_transpose33(P_star_av) /1.e6_pReal math_transpose33(P_av) /1.e6_pReal
lambda_av = lambda_av *wgt lambda_av = lambda_av *wgt
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',&
math_transpose33(lambda_av) /1.e6_pReal math_transpose33(lambda_av) /1.e6_pReal
err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av)) err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av))
err_p = err_p/sqrt(math_mul33xx33(P_star_av,P_star_av)) err_p = err_p/sqrt(math_mul33xx33(P_av,P_av))
write(6,'(a,es14.7,es14.7)') 'error F', err_f/1e-4, err_f write(6,'(a,es14.7,es14.7)') 'error F', err_f/1e-4, err_f
write(6,'(a,es14.7,es14.7)') 'error P', err_p/1e-3, err_p write(6,'(a,es14.7,es14.7)') 'error P', err_p/1e-3, err_p
write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol, err_stress write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol, err_stress
write(6,'(a,es14.7,es14.7)') 'error divergence = ',err_div/err_div_tol, err_div
write(6,*) ' ' write(6,*) ' '
write(6,'(a,es14.7)') 'error divergence FT RMS = ',err_div_RMS
write(6,'(a,es14.7)') 'max abs err F', err_f_point write(6,'(a,es14.7)') 'max abs err F', err_f_point
write(6,'(a,es14.7)') 'max abs err P', err_p_point write(6,'(a,es14.7)') 'max abs err P', err_p_point
err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol) err_crit = max(err_p/1e-3, err_f/1e-4,err_stress/err_stress_tol)
write(6,'(a,es14.7)') 'critical error', err_crit write(6,'(a,es14.7)') 'critical error', err_crit
if (.not. callCPFEM) then
if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true.
err_crit =huge(1.0_pReal)
else
if(iter >2 .and. iter< itmax-3) callCPFEM=.false.
guesses = 0_pInt
endif
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
write(6,'(a)') ' ' write(6,'(a)') ' '
@ -888,15 +876,6 @@ program DAMASK_spectral_AL
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.
write(6,*) 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_star)) ! writing deformation gradient field to file
write (777,rec=1) F_star
close (777)
restartInc=totalIncsCounter
endif
endif ! end calculation/forwarding endif ! end calculation/forwarding
enddo ! end looping over incs in current loadcase enddo ! end looping over incs in current loadcase
deallocate(c_reduced) deallocate(c_reduced)
@ -907,6 +886,6 @@ program DAMASK_spectral_AL
write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', & write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' increments did not converge!' notConvergedCounter + convergedCounter, ' increments did not converge!'
close(538) close(538)
call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) call fftw_destroy_plan(plan_tau); call fftw_destroy_plan(plan_correction)
call quit(1_pInt) call quit(1_pInt)
end program DAMASK_spectral_AL end program DAMASK_spectral_AL

8
code/spectral_quit.f90
View File

@ -53,10 +53,10 @@ subroutine quit(stop_id)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
if (stop_id == 1_pInt) stop 1 ! normal termination if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id <= 0_pInt) then ! trigger regridding if (stop_id < 0_pInt) then ! trigger regridding
write(6,'(a,i6)') 'restart a', stop_id*(-1_pInt) write(0,'(a,i6)') 'restart a', stop_id*(-1_pInt)
stop 2 stop 2
endif endif
stop 0 ! error stop 1 ! error
end subroutine end subroutine