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

added time_sensitive flag to ignore or add time sensitive component to analytic jacobian (should be close to perturbed jacobian when ignored) 

now uses sub_dt for analytic jacobian calculations
bug fixes
This commit is contained in:
Pratheek Shanthraj 2012-03-01 19:42:43 +00:00
parent df03dee91c
commit 17e9698659
3 changed files with 24 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
code/config/numerics.config
View File

@ -11,6 +11,7 @@ integrator 1 # integration method (1 = Fixed Point Iter
integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
Lp_frac 0.5 # when integrating Fp from previous step to current Lp is calculated as (1-Lp_frac)*Lp previous + Lp_frac*Lp current Lp_frac 0.5 # when integrating Fp from previous step to current Lp is calculated as (1-Lp_frac)*Lp previous + Lp_frac*Lp current
analyticJaco 0 # use analytic Jacobian or perturbation (0 = perturbations, 1 = analytic) analyticJaco 0 # use analytic Jacobian or perturbation (0 = perturbations, 1 = analytic)
time_sensitive 0 # adds time sensitive component to analytic jacobian when 1
## crystallite numerical parameters ## ## crystallite numerical parameters ##
nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")

36
code/crystallite.f90
View File

@ -460,7 +460,8 @@ use numerics, only: subStepMinCryst, &
numerics_integrationMode, & numerics_integrationMode, &
relevantStrain, & relevantStrain, &
Lp_frac, & Lp_frac, &
analyticJaco analyticJaco, &
time_sensitive
use debug, only: debug_verbosity, & use debug, only: debug_verbosity, &
debug_selectiveDebugger, & debug_selectiveDebugger, &
debug_e, & debug_e, &
@ -976,14 +977,14 @@ if(updateJaco) then
do g = 1_pInt,myNgrains do g = 1_pInt,myNgrains
Fp_inv_0 = math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)) Fp_inv_0 = math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))
C = math_Mandel66to3333(constitutive_homogenizedC(g,i,e)) C = math_Mandel66to3333(constitutive_homogenizedC(g,i,e))
call math_spectralDecompositionSym33(-(1.0_pReal-Lp_frac)*crystallite_dt(g,i,e)*0.5_pReal& call math_spectralDecompositionSym33(-(1.0_pReal-Lp_frac)*crystallite_subdt(g,i,e)*0.5_pReal&
*crystallite_subLp0(1:3,1:3,g,i,e),Eigval,Eigvec,error) *crystallite_subLp0(1:3,1:3,g,i,e),Eigval,Eigvec,error)
Fp_exp1 = 0.0_pReal Fp_exp1 = 0.0_pReal
Fp_exp1(1,1) = exp(Eigval(1)) Fp_exp1(1,1) = exp(Eigval(1))
Fp_exp1(2,2) = exp(Eigval(2)) Fp_exp1(2,2) = exp(Eigval(2))
Fp_exp1(3,3) = exp(Eigval(3)) Fp_exp1(3,3) = exp(Eigval(3))
Fp_exp1 = math_mul33x33(math_mul33x33(Eigvec,Fp_exp1),math_transpose33(Eigvec)) ! exp(-(1-Lp_frac)*Lp old*dt/2) Fp_exp1 = math_mul33x33(math_mul33x33(Eigvec,Fp_exp1),math_transpose33(Eigvec)) ! exp(-(1-Lp_frac)*Lp old*dt/2)
call math_spectralDecompositionSym33(-(Lp_frac)*crystallite_dt(g,i,e)*0.5_pReal*& call math_spectralDecompositionSym33(-(Lp_frac)*crystallite_subdt(g,i,e)*0.5_pReal*&
crystallite_Lp(1:3,1:3,g,i,e),Eigval,Eigvec,error) crystallite_Lp(1:3,1:3,g,i,e),Eigval,Eigvec,error)
Fp_exp2 = 0.0_pReal Fp_exp2 = 0.0_pReal
Fp_exp2(1,1) = exp(Eigval(1)) Fp_exp2(1,1) = exp(Eigval(1))
@ -1020,30 +1021,26 @@ if(updateJaco) then
Tensor2 = 0.0_pReal Tensor2 = 0.0_pReal
Tensor3 = 0.0_pReal Tensor3 = 0.0_pReal
Tensor4 = 0.0_pReal Tensor4 = 0.0_pReal
Tensor5 = 0.0_pReal
do h=1_pInt,3_pInt; do j=1_pInt,3_pInt do h=1_pInt,3_pInt; do j=1_pInt,3_pInt
V_dir = 0.0_pReal V_dir = 0.0_pReal
V_dir(h,j) = -Lp_frac*crystallite_dt(g,i,e) V_dir(h,j) = -Lp_frac*crystallite_subdt(g,i,e)
V_dir = math_mul33x33(math_mul33x33(math_transpose33(Eigvec),V_dir),Eigvec) V_dir = math_mul33x33(math_mul33x33(math_transpose33(Eigvec),V_dir),Eigvec)
Tensor2(1:3,1:3,h,j) = math_mul33x33(math_mul33x33(dSdFe_mat1,V_dir*phi_mat),dSdFe_mat2) if (time_sensitive) then
Tensor1(h,j,1:3,1:3) = Fpinv_rate(h,j)*FDot_temp ! assuming dF_old = dF_new... need full-rank Tensor1 otherwise. good only for unidirectional loading Tensor1(h,j,1:3,1:3) = Fpinv_rate(h,j)*FDot_temp ! assuming dF_old = dF_new... need full-rank Tensor1 otherwise. good only for unidirectional loading
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt endif
Tensor1(h,j,o,p) = math_I3(h,o)*Fp_inv_current(p,j) + Tensor1(h,j,o,p) ! dF_im/dF_kl * (Fp current^-1)_mj + d(Fp current^-1)_ij/dt * dt/Fdot_kl Tensor1(h,j,h,1:3) = Tensor1(h,j,h,1:3) + Fp_inv_current(1:3,j) ! dF_im/dF_kl * (Fp current^-1)_mj + d(Fp current^-1)_ij/dt * dt/Fdot_kl
Tensor3(h,j,1:3,1:3) = 0.5_pReal*(math_mul33x33(C(h,j,1:3,1:3),math_transpose33(& Tensor2(1:3,1:3,h,j) = math_mul33x33(math_mul33x33(dSdFe_mat1,V_dir*phi_mat),dSdFe_mat2)
crystallite_subFe0(1:3,1:3,g,i,e))) + math_mul33x33(C(h,j,1:3,1:3),& Tensor3(h,j,1:3,1:3) = math_mul33x33(C(h,j,1:3,1:3), &
math_transpose33(crystallite_subFe0(1:3,1:3,g,i,e)))) ! dS_ij/dFe_kl math_transpose33(crystallite_subFe0(1:3,1:3,g,i,e))) ! dS_ij/dFe_kl
Tensor5(h,j,o,p) = math_I3(h,o)*math_I3(j,p) - math_I3(h,j)*math_I3(o,p)/3.0_pReal ! tensor to strip away volumetric components of Lp added due to numerical error enddo; enddo
enddo; enddo; enddo; enddo
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, & call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, &
crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g, i, e) crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
dLp_dT = math_Plain99to3333(dLp_dT_constitutive) dLp_dT = math_Plain99to3333(dLp_dT_constitutive)
Tensor4 = math_mul3333xx3333(math_mul3333xx3333(math_mul3333xx3333(Tensor2,dLp_dT),& Tensor4 = math_mul3333xx3333(math_mul3333xx3333(Tensor2,dLp_dT),Tensor3) ! applying chain rule to get F_im * dFp^-1_mj/dFe_kl
Tensor5),Tensor3) ! applying chain rule to get F_im * dFp^-1_mj/dFe_kl
A99 = math_identity2nd(9_pInt) - math_Plain3333to99(Tensor4) A99 = math_identity2nd(9_pInt) - math_Plain3333to99(Tensor4)
C99 = math_Plain3333to99(Tensor1) ! [I - F_im * dFp^-1_mj/dFe_op]*dFe_op/dF_kl = dF_im/dF_kl * (Fp current^-1)_mj + d(Fp current^-1)_ij/dt * dt/Fdot_kl C99 = math_Plain3333to99(Tensor1) ! [I - F_im * dFp^-1_mj/dFe_op]*dFe_op/dF_kl = dF_im/dF_kl * (Fp current^-1)_mj + d(Fp current^-1)_ij/dt * dt/Fdot_kl
call math_invert(9_pInt,A99, A_inv99, AnzNegEW, error) call math_invert(9_pInt,A99, A_inv99, AnzNegEW, error)
dFedF99 = math_mul99x99(A_inv99,C99) ! solve for dFe_ij/dF_kl dFedF99 = math_mul99x99(A_inv99,C99) ! solve for dFe_ij/dF_kl
dFedF = math_Plain99to3333(dFedF99)
fixedpt_iter = 1_pInt fixedpt_iter = 1_pInt
fixedpt_error = 1.0_pReal fixedpt_error = 1.0_pReal
do while ((fixedpt_iter .lt. 10_pInt) .and. (fixedpt_error .gt. 1e-20_pReal)) ! if solution is not accurate, use as estimate for better solution do while ((fixedpt_iter .lt. 10_pInt) .and. (fixedpt_error .gt. 1e-20_pReal)) ! if solution is not accurate, use as estimate for better solution
@ -1056,18 +1053,19 @@ if(updateJaco) then
fixedpt_iter = fixedpt_iter + 1_pInt fixedpt_iter = fixedpt_iter + 1_pInt
enddo enddo
dFedF = math_Plain99to3333(dFedF99) dFedF = math_Plain99to3333(dFedF99)
Tensor4 = math_mul3333xx3333(Tensor4,dFedF)
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dFedF(1:3,1:3,o,p) = math_transpose33(dFedF(1:3,1:3,o,p)) dFedF(1:3,1:3,o,p) = math_transpose33(dFedF(1:3,1:3,o,p))
enddo; enddo enddo; enddo
Tensor4 = math_mul3333xx3333(Tensor4,dFedF)
dSdF = math_mul3333xx3333(Tensor3,dFedF) ! dS/dF = dS/dFe * dFe/dF dSdF = math_mul3333xx3333(Tensor3,dFedF) ! dS/dF = dS/dFe * dFe/dF
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
crystallite_dPdF(1:3,1:3,o,p,g,i,e) = math_mul33x33(math_mul33x33(math_transpose33(& crystallite_dPdF(1:3,1:3,o,p,g,i,e) = math_mul33x33(math_mul33x33(math_transpose33(&
dFedF(1:3,1:3,o,p)),math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(& dFedF(1:3,1:3,o,p)),math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(&
Fp_inv_current)) + & ! dP/dF = dFe/dF * S * Fp^-T... Fp_inv_current)) + & ! dP/dF = dFe/dF * S * Fp^-T...
math_mul33x33(math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),& math_mul33x33(math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),&
Fp_inv_current),math_Mandel6to33(crystallite_Tstar_v)),math_mul33x33(math_inv33& Fp_inv_current),math_Mandel6to33(crystallite_Tstar_v)),math_transpose33( &
(crystallite_subF(1:3,1:3,g,i,e)),Tensor4(1:3,1:3,o,p))+Fpinv_rate*FDot_temp(o,p))& ! + Fe * S * dFp^-T/dF... math_mul33x33(math_inv33(crystallite_subF(1:3,1:3,g,i,e)),Tensor4(1:3,1:3,o,p) &
+ Fpinv_rate*FDot_temp(o,p)))) & ! + Fe * S * dFp^-T/dF...
+ math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),Fp_inv_current),& + math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),Fp_inv_current),&
math_mul33x33(dSdF(1:3,1:3,o,p),math_transpose33(Fp_inv_current))) ! + Fe * dS/dF * Fp^-T math_mul33x33(dSdF(1:3,1:3,o,p),math_transpose33(Fp_inv_current))) ! + Fe * dS/dF * Fp^-T
enddo; enddo enddo; enddo

8
code/numerics.f90
View File

@ -78,7 +78,8 @@ logical :: memory_efficient = .true., &
divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction
update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness
!* end of spectral parameters: !* end of spectral parameters:
analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations analyticJaco = .false., & ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
time_sensitive = .true. ! adds time sensitive component to analytic jacobian when .true.
@ -200,9 +201,11 @@ subroutine numerics_init()
case ('integratorstiffness') case ('integratorstiffness')
numerics_integrator(2) = IO_intValue(line,positions,2_pInt) numerics_integrator(2) = IO_intValue(line,positions,2_pInt)
case ('lp_frac') case ('lp_frac')
Lp_frac = IO_intValue(line,positions,2_pInt) Lp_frac = IO_floatValue(line,positions,2_pInt)
case ('analyticjaco') case ('analyticjaco')
analyticJaco = IO_intValue(line,positions,2_pInt) > 0_pInt analyticJaco = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('time_sensitive')
time_sensitive = IO_intValue(line,positions,2_pInt) > 0_pInt
!* RGC parameters: !* RGC parameters:
@ -315,6 +318,7 @@ subroutine numerics_init()
write(6,'(a24,2(1x,i8),/)')' integrator: ',numerics_integrator write(6,'(a24,2(1x,i8),/)')' integrator: ',numerics_integrator
write(6,'(a24,1x,e8.1)') ' Lp_frac: ',Lp_frac write(6,'(a24,1x,e8.1)') ' Lp_frac: ',Lp_frac
write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco
write(6,'(a24,1x,L8)') ' time sensitive: ',time_sensitive
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,e8.1)') ' subStepMinHomog: ',subStepMinHomog write(6,'(a24,1x,e8.1)') ' subStepMinHomog: ',subStepMinHomog