From c2f5cebacb02a46f4ebc1b505a778a75dae6c20a Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Wed, 14 Mar 2012 13:56:50 +0000 Subject: [PATCH] simplified analytic jacobian calculation. removed Lpfrac, time_sensitive. introduced rate_sensitivity flag when calling crystallite_stressAndItsTangent that is currently set to .false. and is to be set according to which dPdF the FE solver is asking for --- code/config/numerics.config | 2 - code/crystallite.f90 | 212 ++++++++++++------------------------ code/homogenization.f90 | 4 +- code/numerics.f90 | 14 +-- 4 files changed, 73 insertions(+), 159 deletions(-) diff --git a/code/config/numerics.config b/code/config/numerics.config index 50a09059b..39cfa858f 100644 --- a/code/config/numerics.config +++ b/code/config/numerics.config @@ -9,9 +9,7 @@ pert_Fg 1.0e-7 # deformation gradient perturbation for gr pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central) integrator 1 # integration method (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 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 ## nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") diff --git a/code/crystallite.f90 b/code/crystallite.f90 index ad0418487..9d7d44817 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -386,7 +386,7 @@ crystallite_orientation0 = crystallite_orientation ! Store initial o enddo !$OMP END PARALLEL DO -call crystallite_stressAndItsTangent(.true.) ! request elastic answers +call crystallite_stressAndItsTangent(.true.,.false.) ! request elastic answers crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback ! *** Output to MARC output file *** @@ -452,7 +452,7 @@ end subroutine crystallite_init !******************************************************************** ! calculate stress (P) and tangent (dPdF) for crystallites !******************************************************************** -subroutine crystallite_stressAndItsTangent(updateJaco) +subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) !*** variables and functions from other modules ***! use numerics, only: subStepMinCryst, & @@ -464,9 +464,7 @@ use numerics, only: subStepMinCryst, & numerics_integrator, & numerics_integrationMode, & relevantStrain, & - Lp_frac, & - analyticJaco, & - time_sensitive + analyticJaco use debug, only: debug_what, & debug_crystallite, & debug_levelBasic, & @@ -485,14 +483,8 @@ use math, only: math_inv33, & math_Mandel6to33, & math_Mandel33to6, & math_I3, & - math_Plain3333to99, & - math_Plain99to3333, & - math_mul99x99, & math_Mandel66to3333, & - math_mul3333xx33, & - math_invert, & - math_mul3333xx3333, & - math_spectralDecompositionSym33 + math_mul3333xx3333 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP use mesh, only: mesh_element, & @@ -508,13 +500,11 @@ use constitutive, only: constitutive_sizeState, & constitutive_partionedState0, & constitutive_homogenizedC, & constitutive_dotState, & - constitutive_dotState_backup, & - constitutive_LpAndItsTangent + constitutive_dotState_backup implicit none !*** input variables ***! -logical, intent(in) :: updateJaco ! flag indicating wehther we want to update the Jacobian (stiffness) or not - +logical, intent(in) :: updateJaco, rate_sensitivity ! flag indicating wehther we want to update the Jacobian (stiffness) or not !*** local variables ***! real(pReal) myPert, & ! perturbation with correct sign formerSubStep @@ -552,45 +542,16 @@ integer(pInt) NiterationCrystallite, & ! number of ite logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & convergenceFlag_backup ! local variables used for calculating analytic Jacobian -real(pReal), dimension(3,3):: Fp_exp1, & - Fp_exp2, & - Fp_inv_current, & - dSdFe_mat1, & - dSdFe_mat2, & - Lp_constitutive, & - Eigvec, & - V_dir, & - phi_mat, & - Fpinv_rate, & - FDot_temp, & - Rot_mat, & - Fp0, & - Fp_inv_0, & - Lp0, & - Lp_current, & - Fe0 +real(pReal), dimension(3,3):: Fpinv_rate, & + FDot_inv real(pReal), dimension(3,3,3,3) :: C, & dSdFe, & - dLp_dT, & - Tensor1, & - Tensor2, & - Tensor3, & - Tensor4, & - Tensor5, & dFedF, & - dSdF -real(pReal), dimension(9,9) :: dLp_dT_constitutive, & - A99, & - A_inv99, & - C99, & - dFedF99, & - dFedF_old99 -real(pReal), dimension(3):: Eigval -real(pReal) :: dt, & - fixedpt_error, & - counter -integer(pInt) :: AnzNegEW, & - fixedpt_iter + dFedFdot, & + dSdF, & + dSdFdot, & + dFp_invdFdot +real(pReal) :: counter logical :: error ! --+>> INITIALIZE TO STARTING CONDITION <<+-- @@ -977,109 +938,70 @@ if(updateJaco) then ! --- CALCULATE ANALYTIC dPdF --- - !$OMP PARALLEL DO PRIVATE(myNgrains,Fp_inv_0,C,Eigval,Eigvec,Fp_exp1,Fp_exp2,Fp_inv_current,& - !$OMP dSdFe_mat1,dSdFe_mat2,Fpinv_rate,FDot_temp,counter,phi_mat,Tensor1,Tensor2,Tensor3,Tensor4,& - !$OMP Tensor5,V_dir,dLp_dT_constitutive,dLp_dT,C99,A99,A_inv99,AnzNegEW,error,& - !$OMP dFedF99,dFedF_old99,dFedF,fixedpt_iter,fixedpt_error,dSdF) + !$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,myNgrains,C) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1_pInt,myNgrains - Fp_inv_0 = math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)) - C = math_Mandel66to3333(constitutive_homogenizedC(g,i,e)) - 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) - Fp_exp1 = 0.0_pReal - Fp_exp1(1,1) = exp(Eigval(1)) - Fp_exp1(2,2) = exp(Eigval(2)) - 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) - 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) - Fp_exp2 = 0.0_pReal - Fp_exp2(1,1) = exp(Eigval(1)) - Fp_exp2(2,2) = exp(Eigval(2)) - Fp_exp2(3,3) = exp(Eigval(3)) - Fp_exp2 = math_mul33x33(math_mul33x33(Eigvec,Fp_exp2),math_transpose33(Eigvec)) ! exp(-(Lp_frac)*Lp current*dt/2) - Fp_inv_current = math_mul33x33(math_mul33x33(math_mul33x33(math_mul33x33(Fp_inv_0, & - Fp_exp1),Fp_exp2),Fp_exp2),Fp_exp1) ! Fp current^-1 = Fp old^-1 * exp(-(1-Lp_frac)*Lp old*dt/2) * exp(-(Lp_frac)*Lp current*dt) * exp(-(1-Lp_frac)*Lp old*dt/2) - dSdFe_mat2 = math_mul33x33(math_mul33x33(Fp_exp1,Fp_exp2),Eigvec) - dSdFe_mat1 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),Fp_inv_0),dSdFe_mat2) - dSdFe_mat2 = math_transpose33(dSdFe_mat2) - Fpinv_rate = -math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),math_mul33x33(Fp_inv_current,& - (1.0_pReal-Lp_frac)*crystallite_subLp0(1:3,1:3,g,i,e)+ Lp_frac*crystallite_Lp(1:3,1:3,g,i,e))) ! F * dFp^-1 = F * dFp^-1/dt *dt... dFp may overshoot dF by small ammount as - FDot_temp = crystallite_subF(1:3,1:3,g,i,e) - crystallite_F0(1:3,1:3,g,i,e) ! dF = dFe + dFp is not strictly enforced (can result in small oscillations in dP/dF) + do g = 1_pInt,myNgrains + C = math_Mandel66to3333(constitutive_homogenizedC(g,i,e)) + dFedF = 0.0_pReal + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + dFedF(p,o,o,1:3) = crystallite_invFp(1:3,p,g,i,e) ! dFe_ij/dF_kl = dF_im/dF_kl * (Fp current^-1)_mj + dSdFe(o,p,1:3,1:3) = math_mul33x33(C(o,p,1:3,1:3), & + math_transpose33(crystallite_subFe0(1:3,1:3,g,i,e))) ! dS_ij/dFe_kl + enddo; enddo + dSdF = math_mul3333xx3333(dSdFe,dFedF) ! dS/dF = dS/dFe * dFe/dF + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + crystallite_dPdF(1:3,1:3,o,p,g,i,e) = math_mul33x33(math_mul33x33(dFedF(1:3,1:3,o,p),& + math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(& + crystallite_invFp(1:3,1:3,g,i,e))) & ! dP/dF = dFe/dF * S * Fp^-T... + + math_mul33x33(crystallite_subFe0(1:3,1:3,g,i,e),& + math_mul33x33(dSdF(1:3,1:3,o,p),math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))) ! + Fe * dS/dF * Fp^-T + enddo; enddo + enddo; enddo; enddo + !$OMP END PARALLEL DO + endif + + if (rate_sensitivity) then + !$OMP PARALLEL DO PRIVATE(dFedFdot,dSdFdot,dSdFe,Fpinv_rate,FDot_inv,counter,dFp_invdFdot,C,myNgrains) + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1_pInt,myNgrains + C = math_Mandel66to3333(constitutive_homogenizedC(g,i,e)) + Fpinv_rate = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),crystallite_Lp(1:3,1:3,g,i,e)) ! dFp^-1 = dFp^-1/dt *dt... dFp may overshoot dF by small ammount as + FDot_inv = crystallite_subF(1:3,1:3,g,i,e) - crystallite_F0(1:3,1:3,g,i,e) counter = 0.0_pReal - phi_mat = 0.0_pReal - do h=1_pInt,3_pInt; do j=1_pInt,3_pInt - if (Eigval(h) == Eigval(j)) then - phi_mat(h,j) = 1.0_pReal - else - phi_mat(h,j) = sinh(Eigval(h) - Eigval(j))/(Eigval(h) - Eigval(j)) - endif - if (abs(FDot_temp(h,j)) .lt. relevantStrain) then - counter = counter + 1.0_pReal - FDot_temp(h,j) = 0.0_pReal + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + if (abs(FDot_inv(o,p)) .lt. relevantStrain) then + FDot_inv(o,p) = 0.0_pReal else - FDot_temp(h,j) = crystallite_dt(g,i,e)/FDot_temp(h,j) + counter = counter + 1.0_pReal + FDot_inv(o,p) = crystallite_dt(g,i,e)/FDot_inv(o,p) endif enddo; enddo - if (counter .lt. 9.0_pReal) then - FDot_temp = FDot_temp/(9.0_pReal - counter) - endif - Tensor1 = 0.0_pReal - Tensor2 = 0.0_pReal - Tensor3 = 0.0_pReal - Tensor4 = 0.0_pReal - do h=1_pInt,3_pInt; do j=1_pInt,3_pInt - V_dir = 0.0_pReal - V_dir(h,j) = -Lp_frac*crystallite_subdt(g,i,e) - V_dir = math_mul33x33(math_mul33x33(math_transpose33(Eigvec),V_dir),Eigvec) - 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 - endif - 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 - Tensor2(1:3,1:3,h,j) = math_mul33x33(math_mul33x33(dSdFe_mat1,V_dir*phi_mat),dSdFe_mat2) - 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 + if (counter .gt. 0.0_pReal) FDot_inv = FDot_inv/counter + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + dFp_invdFdot(o,p,1:3,1:3) = Fpinv_rate(o,p)*FDot_inv ! dFe_ij/dF_kl = dF_im/dF_kl * (Fp current^-1)_mj + dSdFe(o,p,1:3,1:3) = math_mul33x33(C(o,p,1:3,1:3), & + math_transpose33(crystallite_subFe0(1:3,1:3,g,i,e))) ! dS_ij/dFe_kl enddo; enddo - call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, & - crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g, i, e) - dLp_dT = math_Plain99to3333(dLp_dT_constitutive) - Tensor4 = math_mul3333xx3333(math_mul3333xx3333(Tensor2,dLp_dT),Tensor3) ! applying chain rule to get F_im * dFp^-1_mj/dFe_kl - 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 - call math_invert(9_pInt,A99, A_inv99, AnzNegEW, error) - dFedF99 = math_mul99x99(A_inv99,C99) ! solve for dFe_ij/dF_kl - fixedpt_iter = 1_pInt - 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 - dFedF_old99 = dFedF99 - A99 = math_mul99x99(A_inv99,A99) - C99 = dFedF_old99 - call math_invert(9_pInt,A99, A_inv99, AnzNegEW, error) - dFedF99 = math_mul99x99(A_inv99,C99) - fixedpt_error = maxval(abs(dFedF99 - dFedF_old99)) - fixedpt_iter = fixedpt_iter + 1_pInt - enddo - dFedF = math_Plain99to3333(dFedF99) - 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)) + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + dFedFdot(1:3,1:3,o,p) = math_transpose33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & + dFp_invdFdot(1:3,1:3,o,p))) enddo; enddo - Tensor4 = math_mul3333xx3333(Tensor4,dFedF) - dSdF = math_mul3333xx3333(Tensor3,dFedF) ! dS/dF = dS/dFe * dFe/dF - 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(& - dFedF(1:3,1:3,o,p)),math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(& - 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),& - Fp_inv_current),math_Mandel6to33(crystallite_Tstar_v)),math_transpose33( & - 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(dSdF(1:3,1:3,o,p),math_transpose33(Fp_inv_current))) ! + Fe * dS/dF * Fp^-T + dSdFdot = math_mul3333xx3333(dSdFe,dFedFdot) + do p=1_pInt,3_pInt; do o=1_pInt,3_pInt + crystallite_dPdF(1:3,1:3,o,p,g,i,e) = crystallite_dPdF(1:3,1:3,o,p,g,i,e) - & + (math_mul33x33(math_mul33x33(dFedFdot(1:3,1:3,o,p), & + math_Mandel6to33(crystallite_Tstar_v)),math_transpose33( & + crystallite_invFp(1:3,1:3,g,i,e))) + & ! dP/dF = dFe/dFdot * S * Fp^-T... + math_mul33x33(math_mul33x33(crystallite_subFe0(1:3,1:3,g,i,e), & + math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(dFp_invdFdot(1:3,1:3,o,p))) & ! + Fe * S * dFp^-T/dFdot... + + math_mul33x33(crystallite_subFe0(1:3,1:3,g,i,e), & + math_mul33x33(dSdFdot(1:3,1:3,o,p),math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))))) ! + Fe * dS/dFdot * Fp^-T enddo; enddo - enddo; enddo; enddo + enddo; enddo; enddo !$OMP END PARALLEL DO endif endif ! jacobian calculation diff --git a/code/homogenization.f90 b/code/homogenization.f90 index ecc6d2806..778cee9fd 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -300,6 +300,7 @@ use debug, only: debug_what, & real(pReal), intent(in) :: dt logical, intent(in) :: updateJaco + logical :: rate_sensitivity integer(pInt) NiterationHomog,NiterationMPstate integer(pInt) g,i,e,myNgrains @@ -485,7 +486,8 @@ use debug, only: debug_what, & ! ! based on crystallite_partionedF0,.._partionedF ! incrementing by crystallite_dt - call crystallite_stressAndItsTangent(updateJaco) ! request stress and tangent calculation for constituent grains + rate_sensitivity = .false. ! request rate sensitive contribution to dPdF + call crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) ! request stress and tangent calculation for constituent grains ! --+>> state update <<+-- diff --git a/code/numerics.f90 b/code/numerics.f90 index 588cef504..f02b33bce 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -51,7 +51,6 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, & rTol_crystalliteTemperature= 1.0e-6_pReal, & ! relative tolerance in crystallite temperature loop rTol_crystalliteStress = 1.0e-6_pReal, & ! relative tolerance in crystallite stress loop aTol_crystalliteStress = 1.0e-8_pReal, & ! absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order - Lp_frac = 0.5_pReal, & ! fraction of Lp current and Lp previous step to use when integrating Fp from previous step to current absTol_RGC = 1.0e+4_pReal, & ! absolute tolerance of RGC residuum relTol_RGC = 1.0e-3_pReal, & ! relative tolerance of RGC residuum @@ -79,8 +78,7 @@ logical :: memory_efficient = .true., & 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 !* end of spectral parameters: - 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. + analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations @@ -195,12 +193,8 @@ subroutine numerics_init numerics_integrator(1) = IO_intValue(line,positions,2_pInt) case ('integratorstiffness') numerics_integrator(2) = IO_intValue(line,positions,2_pInt) - case ('lp_frac') - Lp_frac = IO_floatValue(line,positions,2_pInt) case ('analyticjaco') analyticJaco = IO_intValue(line,positions,2_pInt) > 0_pInt - case ('time_sensitive') - time_sensitive = IO_intValue(line,positions,2_pInt) > 0_pInt !* RGC parameters: @@ -310,10 +304,8 @@ subroutine numerics_init write(6,'(a24,1x,e8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,1x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress - write(6,'(a24,2(1x,i8),/)')' integrator: ',numerics_integrator - write(6,'(a24,1x,e8.1)') ' Lp_frac: ',Lp_frac - write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco - write(6,'(a24,1x,L8)') ' time sensitive: ',time_sensitive + write(6,'(a24,2(1x,i8))')' integrator: ',numerics_integrator + write(6,'(a24,1x,L8,/)') ' analytic Jacobian: ',analyticJaco write(6,'(a24,1x,i8)') ' nHomog: ',nHomog write(6,'(a24,1x,e8.1)') ' subStepMinHomog: ',subStepMinHomog