diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 7abd737f3..399e63b0e 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -518,8 +518,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) nCryst, & numerics_integrator, & numerics_integrationMode, & - numerics_timeSyncing, & - analyticJaco + numerics_timeSyncing use debug, only: & debug_level, & debug_crystallite, & @@ -582,23 +581,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar ! 2nd Piola-Kirchhoff stress tensor - real(pReal), allocatable, dimension(:,:,:,:,:,:,:) :: & - dPdF_perturbation1, & - dPdF_perturbation2 - real(pReal), allocatable, dimension(:,:,:,:,:) :: & - F_backup, & - Fp_backup, & - InvFp_backup, & - Fi_backup, & - InvFi_backup, & - Fe_backup, & - Lp_backup, & - Li_backup, & - P_backup - real(pReal), allocatable, dimension(:,:,:,:) :: & - Tstar_v_backup - logical, allocatable, dimension(:,:,:) :: & - convergenceFlag_backup integer(pInt) :: & NiterationCrystallite, & ! number of iterations in crystallite loop c, & !< counter in integration point component loop @@ -1137,371 +1119,115 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! --+>> STIFFNESS CALCULATION <<+-- computeJacobian: if(updateJaco) then - jacobianMethod: if (analyticJaco) then + !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,& + !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) + elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) + myNcomponents = 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 c = 1_pInt,myNcomponents + call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent - ! --- ANALYTIC JACOBIAN --- + call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e), & + c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration + if (sum(abs(dLidS)) < tol_math_check) then + dFidS = 0.0_pReal + else + temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) + lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal + do o=1_pInt,3_pInt; do p=1_pInt,3_pInt + lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) + & + crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidFi(1:3,1:3,o,p)) + lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) + & + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) + rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) - & + crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidS(1:3,1:3,o,p)) + enddo; enddo + call math_invert(9_pInt,math_Plain3333to99(lhs_3333),temp_99,error) + if (error) then + call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dFidS = 0.0_pReal + else + dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) + endif + dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS + endif - !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,& - !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) - elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = 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 c = 1_pInt,myNcomponents - call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent + call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration + dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS - call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e), & - c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration - if (sum(abs(dLidS)) < tol_math_check) then - dFidS = 0.0_pReal - else - temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) - lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal - do o=1_pInt,3_pInt; do p=1_pInt,3_pInt - lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) + & - crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidFi(1:3,1:3,o,p)) - lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) + & - crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) - rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) - & - crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidS(1:3,1:3,o,p)) - enddo; enddo - call math_invert(9_pInt,math_Plain3333to99(lhs_3333),temp_99,error) - if (error) then - call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dFidS = 0.0_pReal - else - dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) - endif - dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif + temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & + crystallite_invFi(1:3,1:3,c,i,e))) + rhs_3333 = 0.0_pReal + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33) - call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration - dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS + temp_3333 = 0.0_pReal + temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & + math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), & + crystallite_invFi(1:3,1:3,c,i,e)) - temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & - crystallite_invFi(1:3,1:3,c,i,e))) - rhs_3333 = 0.0_pReal - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33) + temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)), & + math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o)) - temp_3333 = 0.0_pReal - temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), & - crystallite_invFi(1:3,1:3,c,i,e)) + lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + & + math_mul3333xx3333(dSdFi,dFidS) - temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o)) + call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error) + if (error) then + call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dSdF = rhs_3333 + else + dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) + endif - lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + & - math_mul3333xx3333(dSdFi,dFidS) + dFpinvdF = 0.0_pReal + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e)* & + math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & + math_mul33x33(temp_3333(1:3,1:3,p,o), & + crystallite_invFi(1:3,1:3,c,i,e))) - call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error) - if (error) then - call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dSdF = rhs_3333 - else - dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) - endif + crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal + temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & + math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), & + math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))) + forall(p=1_pInt:3_pInt) & + crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33) - dFpinvdF = 0.0_pReal - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e)* & - math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & - math_mul33x33(temp_3333(1:3,1:3,p,o), & - crystallite_invFi(1:3,1:3,c,i,e))) + temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), & + math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & + math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33) - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal - temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & - math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), & - math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))) - forall(p=1_pInt:3_pInt) & - crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33) + temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & + math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), & + math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))) - temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), & - math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & - math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33) + temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)), & + math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e))) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & + math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o))) - temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & - math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), & - math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))) - - temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e))) - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & - math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o))) - - enddo; enddo - enddo elementLooping6 - !$OMP END PARALLEL DO - - else jacobianMethod - - ! --- STANDARD (PERTURBATION METHOD) FOR JACOBIAN --- - - numerics_integrationMode = 2_pInt - - ! --- BACKUP --- - allocate(dPdF_perturbation1(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(dPdF_perturbation2(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(F_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Fp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(InvFp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Fi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(InvFi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Fe_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Lp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Li_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(P_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(Tstar_v_backup (6, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal) - allocate(convergenceFlag_backup (homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = .false.) - - !$OMP PARALLEL DO PRIVATE(myNcomponents) - elementLooping7: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents - - plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) - enddo - - plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) - enddo - - F_backup(1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) ! ... and kinematics - Fp_backup(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e) - InvFp_backup(1:3,1:3,c,i,e) = crystallite_invFp(1:3,1:3,c,i,e) - Fi_backup(1:3,1:3,c,i,e) = crystallite_Fi(1:3,1:3,c,i,e) - InvFi_backup(1:3,1:3,c,i,e) = crystallite_invFi(1:3,1:3,c,i,e) - Fe_backup(1:3,1:3,c,i,e) = crystallite_Fe(1:3,1:3,c,i,e) - Lp_backup(1:3,1:3,c,i,e) = crystallite_Lp(1:3,1:3,c,i,e) - Li_backup(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e) - Tstar_v_backup(1:6,c,i,e) = crystallite_Tstar_v(1:6,c,i,e) - P_backup(1:3,1:3,c,i,e) = crystallite_P(1:3,1:3,c,i,e) - convergenceFlag_backup(c,i,e) = crystallite_converged(c,i,e) - enddo; enddo - enddo elementLooping7 - !$END PARALLEL DO - ! --- CALCULATE STATE AND STRESS FOR PERTURBATION --- - - dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment - dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment - pertubationLoop: do perturbation = 1,2 ! forward and backward perturbation - if (iand(pert_method,perturbation) > 0_pInt) then ! mask for desired direction - myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step - do k = 1,3; do l = 1,3 ! ...alter individual components - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & - .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & - .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) & - write(6,'(a,2(1x,i1),1x,a,/)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]' - ! --- INITIALIZE UNPERTURBED STATE --- - - select case(numerics_integrator(numerics_integrationMode)) - case(1_pInt) -!why not OMP? ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents - - plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e)) - enddo - - plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) - enddo - - crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e) - crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e) - crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e) - crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e) - crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e) - crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e) - crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e) - enddo; enddo - enddo - case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step - case(4_pInt,5_pInt) -!why not OMP? ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents - - plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%subState0(:,phasememberAt(c,i,e)) - enddo - - plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) - enddo - - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) - crystallite_Fe(1:3,1:3,c,i,e) = crystallite_subFe0(1:3,1:3,c,i,e) - crystallite_Lp(1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) - crystallite_Li(1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) - crystallite_Tstar_v(1:6,c,i,e) = crystallite_subTstar0_v(1:6,c,i,e) - enddo; enddo - enddo - end select - - ! --- PERTURB EITHER FORWARD OR BACKWARD --- -!why not OMP? - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,myNcomponents - crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e) - crystallite_subF(k,l,c,i,e) = crystallite_subF(k,l,c,i,e) + myPert - crystallite_todo(c,i,e) = crystallite_requested(c,i,e) & - .and. convergenceFlag_backup(c,i,e) - if (crystallite_todo(c,i,e)) crystallite_converged(c,i,e) = .false. ! start out non-converged - enddo; enddo; enddo - - - select case(numerics_integrator(numerics_integrationMode)) - case(1_pInt) - call crystallite_integrateStateFPI() - case(2_pInt) - call crystallite_integrateStateEuler() - case(3_pInt) - call crystallite_integrateStateAdaptiveEuler() - case(4_pInt) - call crystallite_integrateStateRK4() - case(5_pInt) - call crystallite_integrateStateRKCK45() - end select - !why not OMP? - elementLooping8: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - select case(perturbation) - case(1_pInt) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update - dPdF_perturbation1(1:3,1:3,k,l,c,i,e) = & - (crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl - case(2_pInt) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update - dPdF_perturbation2(1:3,1:3,k,l,c,i,e) = & - (crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl - end select - enddo elementLooping8 - - enddo; enddo ! k,l component perturbation loop - - endif - enddo pertubationLoop - - ! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE --- - - elementLooping9: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - select case(pert_method) - case(1_pInt) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 1: central solution converged - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e) - case(2_pInt) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 2: central solution converged - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e) - case(3_pInt) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 3: central solution converged - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e) & - + dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e)) - end select - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, & - crystallite_requested(c,i,e) .and. .not. convergenceFlag_backup(c,i,e)) & ! for any pertubation mode: if central solution did not converge... - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,c,i,e) ! ...use (elastic) fallback - enddo elementLooping9 - - ! --- RESTORE --- -!why not OMP? - elementLooping10: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents - - plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e)) - enddo - - plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) - enddo - - crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e) - crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e) - crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e) - crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e) - crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e) - crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e) - crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e) - crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e) - crystallite_P(1:3,1:3,c,i,e) = P_backup(1:3,1:3,c,i,e) - crystallite_converged(c,i,e) = convergenceFlag_backup(c,i,e) - enddo; enddo - enddo elementLooping10 - - deallocate(dPdF_perturbation1) - deallocate(dPdF_perturbation2) - deallocate(F_backup ) - deallocate(Fp_backup ) - deallocate(InvFp_backup ) - deallocate(Fi_backup ) - deallocate(InvFi_backup ) - deallocate(Fe_backup ) - deallocate(Lp_backup ) - deallocate(Li_backup ) - deallocate(P_backup ) - deallocate(Tstar_v_backup ) - deallocate(convergenceFlag_backup) - - endif jacobianMethod + enddo; enddo + enddo elementLooping6 + !$OMP END PARALLEL DO endif computeJacobian !why not OMP? diff --git a/code/numerics.f90 b/code/numerics.f90 index 9b3449e95..9840630db 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -64,7 +64,6 @@ module numerics charLength = 1.0_pReal, & !< characteristic length scale for gradient problems residualStiffness = 1.0e-6_pReal !< non-zero residual damage logical, protected, public :: & - analyticJaco = .true., & !< use analytic Jacobian or perturbation, Default for Spectral solver .true.: usePingPong = .true., & numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity @@ -315,8 +314,6 @@ subroutine numerics_init numerics_integrator(1) = IO_intValue(line,chunkPos,2_pInt) case ('integratorstiffness') numerics_integrator(2) = IO_intValue(line,chunkPos,2_pInt) - case ('analyticjaco') - analyticJaco = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('usepingpong') usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('timesyncing') @@ -528,7 +525,6 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing - write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 index 0e7105c70..e91883173 100644 --- a/code/plastic_disloUCLA.f90 +++ b/code/plastic_disloUCLA.f90 @@ -152,8 +152,6 @@ subroutine plastic_disloUCLA_init(fileUnit) MATERIAL_partPhase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -173,11 +171,9 @@ subroutine plastic_disloUCLA_init(fileUnit) line = '' real(pReal), dimension(:), allocatable :: tempPerSlip - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt) if (maxNinstance == 0_pInt) return @@ -498,10 +494,6 @@ subroutine plastic_disloUCLA_init(fileUnit) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 0902ecad6..e205f5abb 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -238,8 +238,6 @@ subroutine plastic_dislotwin_init(fileUnit) MATERIAL_partPhase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -261,11 +259,9 @@ subroutine plastic_dislotwin_init(fileUnit) line = '' real(pReal), dimension(:), allocatable :: tempPerSlip, tempPerTwin, tempPerTrans - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOTWIN_ID),pInt) if (maxNinstance == 0_pInt) return @@ -930,10 +926,6 @@ subroutine plastic_dislotwin_init(fileUnit) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_isotropic.f90 b/code/plastic_isotropic.f90 index 7d3e1aa7c..9f8129fac 100644 --- a/code/plastic_isotropic.f90 +++ b/code/plastic_isotropic.f90 @@ -94,8 +94,6 @@ subroutine plastic_isotropic_init(fileUnit) debug_constitutive, & debug_levelBasic use numerics, only: & - analyticJaco, & - worldrank, & numerics_integrator use math, only: & math_Mandel3333to66, & @@ -145,11 +143,9 @@ subroutine plastic_isotropic_init(fileUnit) outputtag = '' integer(pInt) :: NipcMyPhase - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt) if (maxNinstance == 0_pInt) return @@ -316,10 +312,6 @@ subroutine plastic_isotropic_init(fileUnit) allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/code/plastic_none.f90 b/code/plastic_none.f90 index a9007667f..7a7589774 100644 --- a/code/plastic_none.f90 +++ b/code/plastic_none.f90 @@ -34,7 +34,6 @@ subroutine plastic_none_init use IO, only: & IO_timeStamp use numerics, only: & - worldrank, & numerics_integrator use material, only: & phase_plasticity, & @@ -53,11 +52,9 @@ subroutine plastic_none_init sizeDotState, & sizeDeltaState - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt) if (maxNinstance == 0_pInt) return @@ -84,11 +81,9 @@ subroutine plastic_none_init allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase)) allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase)) allocate(plasticState(phase)%state (sizeState,NofMyPhase)) - allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase)) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase)) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase)) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase)) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase)) allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase)) diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index a6c7acad1..cb2b31772 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -295,8 +295,6 @@ use material, only: phase_plasticity, & material_phase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator @@ -332,11 +330,9 @@ integer(pInt) :: phase, & integer(pInt) :: NofMyPhase - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) if (maxNinstances == 0) return ! we don't have to do anything if there's no instance for this constitutive law @@ -1310,10 +1306,6 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_phenoplus.f90 b/code/plastic_phenoplus.f90 index 5df3d490b..b9d09fb3d 100644 --- a/code/plastic_phenoplus.f90 +++ b/code/plastic_phenoplus.f90 @@ -145,8 +145,6 @@ subroutine plastic_phenoplus_init(fileUnit) MATERIAL_partPhase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -168,11 +166,9 @@ subroutine plastic_phenoplus_init(fileUnit) line = '' real(pReal), dimension(:), allocatable :: tempPerSlip - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPLUS_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPLUS_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPLUS_ID),pInt) if (maxNinstance == 0_pInt) return @@ -589,10 +585,6 @@ subroutine plastic_phenoplus_init(fileUnit) allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/code/plastic_phenopowerlaw.f90 b/code/plastic_phenopowerlaw.f90 index 12debbb23..0db6a6072 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -157,8 +157,6 @@ subroutine plastic_phenopowerlaw_init(fileUnit) MATERIAL_partPhase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -181,11 +179,9 @@ subroutine plastic_phenopowerlaw_init(fileUnit) line = '' real(pReal), dimension(:), allocatable :: tempPerSlip - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt) if (maxNinstance == 0_pInt) return @@ -587,10 +583,6 @@ subroutine plastic_phenopowerlaw_init(fileUnit) allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/code/plastic_titanmod.f90 b/code/plastic_titanmod.f90 index 24bf543b7..ac80af82b 100644 --- a/code/plastic_titanmod.f90 +++ b/code/plastic_titanmod.f90 @@ -216,8 +216,6 @@ subroutine plastic_titanmod_init(fileUnit) MATERIAL_partPhase use lattice use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -241,11 +239,9 @@ subroutine plastic_titanmod_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_TITANMOD_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_TITANMOD_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_plasticity == PLASTICITY_TITANMOD_ID),pInt) if (maxNinstance == 0_pInt) return @@ -859,10 +855,6 @@ subroutine plastic_titanmod_init(fileUnit) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/prec.f90 b/code/prec.f90 index 0ded23528..190aca02b 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -67,11 +67,9 @@ module prec real(pReal), allocatable, dimension(:,:) :: & partionedState0, & subState0, & - state_backup, & deltaState, & previousDotState, & !< state rate of previous xxxx previousDotState2, & !< state rate two xxxx ago - dotState_backup, & !< backup of state rate RK4dotState real(pReal), allocatable, dimension(:,:,:) :: & RKCK45dotState diff --git a/code/source_damage_anisoBrittle.f90 b/code/source_damage_anisoBrittle.f90 index efd76091e..53cc411af 100644 --- a/code/source_damage_anisoBrittle.f90 +++ b/code/source_damage_anisoBrittle.f90 @@ -92,8 +92,6 @@ subroutine source_damage_anisoBrittle_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator use lattice, only: & lattice_maxNcleavageFamily, & @@ -111,11 +109,9 @@ subroutine source_damage_anisoBrittle_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoBrittle_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoBrittle_LABEL//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt) if (maxNinstance == 0_pInt) return @@ -268,10 +264,6 @@ subroutine source_damage_anisoBrittle_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_anisoDuctile.f90 b/code/source_damage_anisoDuctile.f90 index 72480382a..1a79e3b34 100644 --- a/code/source_damage_anisoDuctile.f90 +++ b/code/source_damage_anisoDuctile.f90 @@ -96,8 +96,6 @@ subroutine source_damage_anisoDuctile_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator use lattice, only: & lattice_maxNslipFamily, & @@ -115,11 +113,9 @@ subroutine source_damage_anisoDuctile_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoDuctile_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoDuctile_LABEL//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt) if (maxNinstance == 0_pInt) return @@ -270,10 +266,6 @@ subroutine source_damage_anisoDuctile_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_isoBrittle.f90 b/code/source_damage_isoBrittle.f90 index 1603ecf48..18194618e 100644 --- a/code/source_damage_isoBrittle.f90 +++ b/code/source_damage_isoBrittle.f90 @@ -82,8 +82,6 @@ subroutine source_damage_isoBrittle_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -97,11 +95,9 @@ subroutine source_damage_isoBrittle_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoBrittle_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoBrittle_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt) if (maxNinstance == 0_pInt) return @@ -222,10 +218,6 @@ subroutine source_damage_isoBrittle_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_isoDuctile.f90 b/code/source_damage_isoDuctile.f90 index 3a85bab24..f30f9a72e 100644 --- a/code/source_damage_isoDuctile.f90 +++ b/code/source_damage_isoDuctile.f90 @@ -82,8 +82,6 @@ subroutine source_damage_isoDuctile_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -97,11 +95,9 @@ subroutine source_damage_isoDuctile_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoDuctile_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoDuctile_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt) if (maxNinstance == 0_pInt) return @@ -222,10 +218,6 @@ subroutine source_damage_isoDuctile_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_thermal_dissipation.f90 b/code/source_thermal_dissipation.f90 index dd6453db0..d649549ad 100644 --- a/code/source_thermal_dissipation.f90 +++ b/code/source_thermal_dissipation.f90 @@ -68,8 +68,6 @@ subroutine source_thermal_dissipation_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -83,11 +81,9 @@ subroutine source_thermal_dissipation_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) if (maxNinstance == 0_pInt) return @@ -163,10 +159,6 @@ subroutine source_thermal_dissipation_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_thermal_externalheat.f90 b/code/source_thermal_externalheat.f90 index 203826205..e2b14d45d 100644 --- a/code/source_thermal_externalheat.f90 +++ b/code/source_thermal_externalheat.f90 @@ -73,8 +73,6 @@ subroutine source_thermal_externalheat_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -89,11 +87,9 @@ subroutine source_thermal_externalheat_init(fileUnit) line = '' real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt) if (maxNinstance == 0_pInt) return @@ -189,10 +185,6 @@ subroutine source_thermal_externalheat_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_irradiation.f90 b/code/source_vacancy_irradiation.f90 index fd7220020..986c229ff 100644 --- a/code/source_vacancy_irradiation.f90 +++ b/code/source_vacancy_irradiation.f90 @@ -70,8 +70,6 @@ subroutine source_vacancy_irradiation_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -85,11 +83,9 @@ subroutine source_vacancy_irradiation_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt) if (maxNinstance == 0_pInt) return @@ -169,10 +165,6 @@ subroutine source_vacancy_irradiation_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_phenoplasticity.f90 b/code/source_vacancy_phenoplasticity.f90 index 2690cf691..924490637 100644 --- a/code/source_vacancy_phenoplasticity.f90 +++ b/code/source_vacancy_phenoplasticity.f90 @@ -68,8 +68,6 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -83,11 +81,9 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt) if (maxNinstance == 0_pInt) return @@ -163,10 +159,6 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_thermalfluc.f90 b/code/source_vacancy_thermalfluc.f90 index 5891ff764..b835e8bce 100644 --- a/code/source_vacancy_thermalfluc.f90 +++ b/code/source_vacancy_thermalfluc.f90 @@ -72,8 +72,6 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) sourceState, & MATERIAL_partPhase use numerics,only: & - analyticJaco, & - worldrank, & numerics_integrator implicit none @@ -87,11 +85,9 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt) if (maxNinstance == 0_pInt) return @@ -170,10 +166,6 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (.not. analyticJaco) then - allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) - endif if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)