From 1195593233c565c39ba61af79bfb1ea30686ef2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 19:39:09 +0200 Subject: [PATCH 01/15] some casting was needed, corrected header --- processing/pre/geom_clean.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/processing/pre/geom_clean.py b/processing/pre/geom_clean.py index c61cf5ef9..408178c32 100755 --- a/processing/pre/geom_clean.py +++ b/processing/pre/geom_clean.py @@ -11,7 +11,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def mostFrequent(arr): - return np.argmax(np.bincount(arr)) + return np.argmax(np.bincount(arr.astype('int'))) #-------------------------------------------------------------------------------------------------- @@ -72,7 +72,13 @@ for name in filenames: # --- do work ------------------------------------------------------------------------------------ - microstructure = ndimage.filters.generic_filter(microstructure,mostFrequent,size=(options.stencil,)*3) + microstructure = ndimage.filters.generic_filter(microstructure,mostFrequent,size=(options.stencil,)*3).astype('int_') + newInfo = {'microstructures': microstructure.max()} + +# --- report --------------------------------------------------------------------------------------- + if ( newInfo['microstructures'] != info['microstructures']): + damask.util.croak('--> microstructures: %i'%newInfo['microstructures']) + info['microstructures'] == newInfo['microstructures'] # --- write header --------------------------------------------------------------------------------- From e4c590699f6ca9743b47e131e4a4f0bcb27b2799 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 20:07:12 +0200 Subject: [PATCH 02/15] removed numerical perturbation calculation --- code/crystallite.f90 | 468 +++++++++---------------------------------- 1 file changed, 97 insertions(+), 371 deletions(-) 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? From faf44c07540739357a264b49f601df743e16051b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 20:12:00 +0200 Subject: [PATCH 03/15] removed allocation for (dot)state backup, simplified screen output handling for MPI --- code/plastic_disloUCLA.f90 | 12 ++---------- code/plastic_dislotwin.f90 | 12 ++---------- code/plastic_isotropic.f90 | 12 ++---------- code/plastic_none.f90 | 7 ++----- code/plastic_nonlocal.f90 | 12 ++---------- code/plastic_phenoplus.f90 | 12 ++---------- code/plastic_phenopowerlaw.f90 | 12 ++---------- code/plastic_titanmod.f90 | 12 ++---------- code/source_damage_anisoBrittle.f90 | 12 ++---------- code/source_damage_anisoDuctile.f90 | 12 ++---------- code/source_damage_isoBrittle.f90 | 12 ++---------- code/source_damage_isoDuctile.f90 | 12 ++---------- code/source_thermal_dissipation.f90 | 12 ++---------- code/source_thermal_externalheat.f90 | 12 ++---------- code/source_vacancy_irradiation.f90 | 12 ++---------- code/source_vacancy_phenoplasticity.f90 | 12 ++---------- code/source_vacancy_thermalfluc.f90 | 12 ++---------- 17 files changed, 34 insertions(+), 165 deletions(-) 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..34570ceab 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 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/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) From d1b777a47c4c54f92260d20a97620945a957610c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 20:12:56 +0200 Subject: [PATCH 04/15] flag for analytic tangent not needed anymore --- code/numerics.f90 | 4 ---- 1 file changed, 4 deletions(-) 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 From d0544a69dd88b5bce04c334ea30b83a27ac5b627 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 20:16:27 +0200 Subject: [PATCH 05/15] removed some remainders of state perturbation --- code/plastic_none.f90 | 2 -- code/prec.f90 | 2 -- 2 files changed, 4 deletions(-) diff --git a/code/plastic_none.f90 b/code/plastic_none.f90 index 34570ceab..7a7589774 100644 --- a/code/plastic_none.f90 +++ b/code/plastic_none.f90 @@ -81,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/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 From 23541a713db1844c839d4a4f1345bbb3140f2f73 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jul 2016 22:26:18 +0200 Subject: [PATCH 06/15] material parameters from paper --- .../ConfigFiles/Phase_Nonlocal_Nickel.config | 146 ++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100755 examples/ConfigFiles/Phase_Nonlocal_Nickel.config diff --git a/examples/ConfigFiles/Phase_Nonlocal_Nickel.config b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config new file mode 100755 index 000000000..1e9ce28fa --- /dev/null +++ b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config @@ -0,0 +1,146 @@ +##################### +# $Id: material.config 477 2009-12-15 10:22:17Z MPIE\c.kords $ +##################### + +##################### + +[SX] +type isostrain +Ngrains 1 + + + +##################### + +##################### + +[Cu_nonlocal] +crystallite 1 +(constituent) phase 1 texture 1 fraction 1.0 + + +##################### + +##################### + +[ori] +(output) eulerangles +(output) grainrotation +(output) grainrotationX +(output) grainrotationY +(output) grainrotationZ +(output) f +(output) fe +(output) p +(output) ee # elastic strain as Green-Lagrange tensor + + +##################### + +##################### + +/echo/ + +[Ni_nonlocal] + +elasticity hooke +plasticity nonlocal +/nonlocal/ + +(output) rho +(output) rho_sgl_mobile +(output) rho_sgl_immobile +(output) rho_sgl_edge_pos +(output) rho_sgl_edge_neg +(output) rho_sgl_screw_pos +(output) rho_sgl_screw_neg +(output) rho_dip_edge +(output) rho_dip_screw +(output) rho_forest +(output) excess_rho_edge +(output) excess_rho_screw +(output) accumulatedshear +(output) shearrate +(output) resolvedstress +(output) resistance +(output) velocity_edge_pos +(output) rho_dot_gen +(output) rho_dot_sgl2dip_edge +(output) rho_dot_sgl2dip_screw +(output) rho_dot_ann_ath +(output) rho_dot_ann_the_edge +(output) rho_dot_ann_the_screw +(output) rho_dot_edgejogs +(output) rho_dot_flux_edge +(output) rho_dot_flux_screw +(output) slipdirection.x +(output) slipdirection.y +(output) slipdirection.z +(output) slipnormal.x +(output) slipnormal.y +(output) slipnormal.z +(output) fluxdensity_edge_pos.x +(output) fluxdensity_edge_pos.y +(output) fluxdensity_edge_pos.z +(output) fluxdensity_edge_neg.x +(output) fluxdensity_edge_neg.y +(output) fluxdensity_edge_neg.z +(output) fluxdensity_screw_pos.x +(output) fluxdensity_screw_pos.y +(output) fluxdensity_screw_pos.z +(output) fluxdensity_screw_neg.x +(output) fluxdensity_screw_neg.y +(output) fluxdensity_screw_neg.z + +lattice_structure fcc +Nslip 12 # number of slip systems per family +c11 246.5e9 +c12 147.3e9 +c44 124.7e9 +burgers 2.48e-10 0 0 0 # Burgers vector in m +rhoSglEdgePos0 6e10 # Initial positive edge single dislocation density in m/m**3 +rhoSglEdgeNeg0 6e10 # Initial negative edge single dislocation density in m/m**3 +rhoSglScrewPos0 6e10 # Initial positive screw single dislocation density in m/m**3 +rhoSglScrewNeg0 6e10 # Initial negative screw single dislocation density in m/m**3 +rhoDipEdge0 0 # Initial edge dipole dislocation density in m/m**3 +rhoDipScrew0 0 # Initial screw dipole dislocation density in m/m**3 +rhoSglScatter 0 +minimumDipoleHeightEdge 2.6e-9 # 3.0e-9 # minimum distance for stable edge dipoles in m +minimumDipoleHeightScrew 12.0e-9 # 50e-9 # minimum distance for stable screw dipoles in m +lambda0 45 # 33 # prefactor for mean free path +edgeMultiplication 0.1 +randomMultiplication 0 +atomicVolume 1.2e-29 +selfdiffusionPrefactor 1.9e-4 # Gottstein p.168 # prefactor for self-diffusion coefficient +selfdiffusionEnergy 5.1e-19 # Gottstein p.168 # activation energy self-diffusion +solidSolutionEnergy 1.8e-19 # activation energy of solid solution particles in J +solidSolutionConcentration 5e-7 # 1e-7 +solidSolutionSize 1.0 +peierlsStressEdge 1e5 # Peierls stress for edges in Pa (per slip family) +peierlsStressScrew 1e5 # Peierls stress for screws in Pa (per slip family) +doublekinkWidth 10 # width of double kinks in multiples of burgers vector length b +viscosity 1e-3 # viscosity for dislocation glide in Pa s +p 1 # exponent for thermal barrier profile +q 1 # exponent for thermal barrier profile +attackFrequency 50e9 # attack frequency in Hz +surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux +grainBoundaryTransmissivity 0.0 +aTol_rho 1e100 # absolute tolerance for dislocation density in m/m**3 +aTol_shear 1e10 # absolute tolerance for dislocation density in m/m**3 +significantRho 1e8 # dislocation density considered relevant in m/m**3 +significantN 1 +shortRangeStressCorrection 0 +CFLfactor 1.1 # safety factor for CFL flux check (numerical parameter) +r 1 +interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient +linetension 0.8 +edgejog 0.01 # 0.2 + + +##################### + +##################### + +[111] +(gauss) phi1 315 Phi 0 phi2 0 scatter 0.000 fraction 1.000 + From 10897bcaaeac3fad132a46457f1941b50d5615bb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 Jul 2016 09:18:13 +0200 Subject: [PATCH 07/15] not needed anymore --- code/damask.core.pyf | 126 ------------------------------------------- 1 file changed, 126 deletions(-) delete mode 100644 code/damask.core.pyf diff --git a/code/damask.core.pyf b/code/damask.core.pyf deleted file mode 100644 index e6396ee1d..000000000 --- a/code/damask.core.pyf +++ /dev/null @@ -1,126 +0,0 @@ -! $Id$ -! -*- f90 -*- -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Note: the syntax of this file is case sensitive. -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! This file was auto-generated with f2py (version:2_5972). -! See http://cens.ioc.ee/projects/f2py2e/ -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! The auto-generated file is quite heavily corrected -! For modifying, notice the following hints: -! - if the dimension of an array depend on a array that is itself an input, use the C-Syntax: (1) becomes [0] etc. -! - be sure that the precision defined is integer, real*8, and complex*16 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -python module core ! in - interface ! in :core - - module prec - subroutine prec_init - end subroutine prec_init - end module prec - - module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90 - subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90 - character(len=1024), intent(in) :: loadcaseParameterIn - character(len=1024), intent(in) :: geometryParameterIn - end subroutine DAMASK_interface_init - end module damask_interface - - module io - subroutine IO_init - end subroutine IO_init - end module io - - module numerics - subroutine numerics_init - end subroutine numerics_init - end module numerics - - module debug - subroutine debug_init - end subroutine debug_init - end module debug - - module math ! in :math:math.f90 - subroutine math_init - end subroutine math_init - - function math_tensorAvg(field) ! in :math:math.f90 - ! input variables - real*8 dimension(:,:,:,:,:), intent(in), :: field - ! function definition - real*8 dimension(3,3), :: math_tensorAvg - end function math_tensorAvg - - end module math - - module fesolving - subroutine FE_init - end subroutine FE_init - end module fesolving - - module mesh ! in :mesh:mesh.f90 - subroutine mesh_init(ip,element) - integer, parameter :: ip = 1 - integer, parameter :: element = 1 - end subroutine mesh_init - - function mesh_nodesAroundCentres(gDim,Favg,centres) ! in :mesh:mesh.f90 - real*8, dimension(:,:,:,:), intent(in) :: centres - real*8, dimension(3), intent(in) :: gDim - real*8, dimension(3,3), intent(in) :: Favg - real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: mesh_nodesAroundCentres - real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: wrappedCentres - end function mesh_nodesAroundCentres - - function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) ! in :mesh:mesh.f90 - real*8, dimension(:,:,:,:,:), intent(in) :: F - real*8, dimension(3), intent(in) :: gDim - real*8, dimension(3,3), intent(in), optional :: FavgIn = -1.0 - real*8, dimension(3), intent(in), optional :: scalingIn = -1.0 - real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsFFT - end function mesh_deformedCoordsFFT - - function mesh_volumeMismatch(gDim,F,nodes) ! in :mesh:mesh.f90 - real*8, dimension(:,:,:,:,:), intent(in) :: F - real*8, dimension(:,:,:,:), intent(in) :: nodes - real*8, dimension(3), intent(in) :: gDim - real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_volumeMismatch - end function mesh_volumeMismatch - - function mesh_shapeMismatch(gDim,F,nodes,centres) ! in :mesh:mesh.f90 - real*8, dimension(:,:,:,:,:), intent(in) :: F - real*8, dimension(:,:,:,:), intent(in) :: nodes - real*8, dimension(:,:,:,:), intent(in) :: centres - real*8, dimension(3), intent(in) :: gDim - real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_shapeMismatch - end function mesh_shapeMismatch - - function mesh_init_postprocessing(filepath) ! in :mesh:mesh.f90 - character(len=*), intent(in) :: filepath - end function mesh_init_postprocessing - - function mesh_build_cellnodes(nodes,Ncellnodes) ! in :mesh:mesh.f90 - integer, intent(in) :: Ncellnodes - real*8, dimension(3,:), intent(in) :: nodes - real*8, dimension(3,Ncellnodes), depend(Ncellnodes) :: mesh_build_cellnodes - end function mesh_build_cellnodes - - function mesh_get_Ncellnodes() ! in :mesh:mesh.f90 - integer :: mesh_get_Ncellnodes - end function mesh_get_Ncellnodes - - function mesh_get_unitlength() ! in :mesh:mesh.f90 - real*8 :: mesh_get_unitlength - end function mesh_get_unitlength - - function mesh_get_nodeAtIP(elemtypeFE,ip) ! in :mesh:mesh.f90 - character(len=*), intent(in) :: elemtypeFE - integer, intent(in) :: ip - integer :: mesh_get_nodeAtIP - end function mesh_get_nodeAtIP - - end module mesh - end interface -end python module core - From 920cf2c849a1e5a1022f72fe3a51e4f834aa04d1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 Jul 2016 09:20:06 +0200 Subject: [PATCH 08/15] removed non-phase parts --- .../ConfigFiles/Phase_Nonlocal_Nickel.config | 52 ------------------- 1 file changed, 52 deletions(-) diff --git a/examples/ConfigFiles/Phase_Nonlocal_Nickel.config b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config index 1e9ce28fa..4d4c2d1df 100755 --- a/examples/ConfigFiles/Phase_Nonlocal_Nickel.config +++ b/examples/ConfigFiles/Phase_Nonlocal_Nickel.config @@ -1,46 +1,3 @@ -##################### -# $Id: material.config 477 2009-12-15 10:22:17Z MPIE\c.kords $ -##################### - -##################### - -[SX] -type isostrain -Ngrains 1 - - - -##################### - -##################### - -[Cu_nonlocal] -crystallite 1 -(constituent) phase 1 texture 1 fraction 1.0 - - -##################### - -##################### - -[ori] -(output) eulerangles -(output) grainrotation -(output) grainrotationX -(output) grainrotationY -(output) grainrotationZ -(output) f -(output) fe -(output) p -(output) ee # elastic strain as Green-Lagrange tensor - - -##################### - -##################### - -/echo/ - [Ni_nonlocal] elasticity hooke @@ -135,12 +92,3 @@ r 1 interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient linetension 0.8 edgejog 0.01 # 0.2 - - -##################### - -##################### - -[111] -(gauss) phi1 315 Phi 0 phi2 0 scatter 0.000 fraction 1.000 - From 360daf9394f7bdc891be664979d9aa218d1ed4d7 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 28 Jul 2016 16:43:38 +0200 Subject: [PATCH 09/15] updated version information after successful test of v2.0.1-5-g920cf2c --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 0ac852dde..5af1fb958 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1 +v2.0.1-5-g920cf2c From 0aadda4cc3c46db849a9e1e94bebe2558f02548b Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 29 Jul 2016 17:45:37 -0400 Subject: [PATCH 10/15] improved cross-platform robustness of DAMASK_env(s) --- DAMASK_env.csh | 28 ++++++++++++++++------------ DAMASK_env.sh | 46 +++++++++++++++++++++++++--------------------- DAMASK_env.zsh | 46 ++++++++++++++++++++-------------------------- 3 files changed, 61 insertions(+), 59 deletions(-) diff --git a/DAMASK_env.csh b/DAMASK_env.csh index 7e3bf128d..78f8922cd 100644 --- a/DAMASK_env.csh +++ b/DAMASK_env.csh @@ -23,25 +23,29 @@ endif # according to http://software.intel.com/en-us/forums/topic/501500 # this seems to make sense for the stack size -set FREE=`which free` -set FREE='' -if ( "x$FREE" != "x" ) then +if ( `which free` != "free: Command not found." ) then set freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'` - set heap=`expr $freeMem / 2` set stack=`expr $freeMem / $DAMASK_NUM_THREADS / 2` + set heap=` expr $freeMem / 2` # http://superuser.com/questions/220059/what-parameters-has-ulimit limit stacksize $stack # maximum stack size (kB) - limit datasize $heap # maximum heap size (kB) + limit datasize $heap # maximum heap size (kB) +endif +if ( `limit | grep coredumpsize` != "" ) then + limit coredumpsize 0 # prevent core dumping +endif +if ( `limit | grep memoryuse` != "" ) then + limit memoryuse unlimited # maximum physical memory size +endif +if ( `limit | grep vmemoryuse` != "" ) then + limit vmemoryuse unlimited # maximum virtual memory size endif -limit coredumpsize 2000 # core file size (512-byte blocks) -limit vmemoryuse unlimited # maximum virtual memory size -limit memoryuse unlimited # maximum physical memory size # disable output in case of scp -if($?prompt) then +if ( $?prompt ) then echo '' - echo Düsseldorf Advanced Materials Simulation Kit - DAMASK - echo Max-Planck-Institut für Eisenforschung, Düsseldorf + echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK + echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf echo https://damask.mpie.de echo echo Using environment with ... @@ -59,8 +63,8 @@ if($?prompt) then echo "MSC.Marc/Mentat $MSC_ROOT" endif echo - echo `limit stacksize` echo `limit datasize` + echo `limit stacksize` endif setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS diff --git a/DAMASK_env.sh b/DAMASK_env.sh index 8536a36ca..d9b5b004f 100644 --- a/DAMASK_env.sh +++ b/DAMASK_env.sh @@ -5,8 +5,8 @@ if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then DAMASK_ROOT=$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$(dirname $BASH_SOURCE)") else - [[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="`pwd`/" - STAT=$(stat "`dirname $BASE$BASH_SOURCE`") + [[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="$(pwd)/" + STAT=$(stat "$(dirname $BASE$BASH_SOURCE)") DAMASK_ROOT=${STAT##* } fi @@ -17,16 +17,16 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set -# if DAMASK_BIN is present and not in $PATH, add it -if [[ "x$DAMASK_BIN" != "x" && ! `echo ":$PATH:" | grep $DAMASK_BIN:` ]]; then +# add DAMASK_BIN if present but not yet in $PATH +if [[ "x$DAMASK_BIN" != "x" && ! $(echo ":$PATH:" | grep $DAMASK_BIN:) ]]; then export PATH=$DAMASK_BIN:$PATH fi -SOLVER=`which DAMASK_spectral 2>/dev/null` +SOLVER=$(which DAMASK_spectral 2>/dev/null) if [ "x$SOLVER" == "x" ]; then SOLVER='Not found!' fi -PROCESSING=`which postResults 2>/dev/null` +PROCESSING=$(which postResults 2>/dev/null) if [ "x$PROCESSING" == "x" ]; then PROCESSING='Not found!' fi @@ -36,25 +36,22 @@ fi # according to http://software.intel.com/en-us/forums/topic/501500 # this seems to make sense for the stack size -FREE=`which free 2>/dev/null` +FREE=$(which free 2>/dev/null) if [ "x$FREE" != "x" ]; then - freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'` - heap=`expr $freeMem / 2` - stack=`expr $freeMem / $DAMASK_NUM_THREADS / 2` - + freeMem=$(free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}') # http://superuser.com/questions/220059/what-parameters-has-ulimit - ulimit -s $stack 2>/dev/null # maximum stack size (kB) - ulimit -d $heap 2>/dev/null # maximum heap size (kB) + ulimit -s $(expr $freeMem / $DAMASK_NUM_THREADS / 2) 2>/dev/null # maximum stack size (kB) + ulimit -d $(expr $freeMem / 2) 2>/dev/null # maximum heap size (kB) fi -ulimit -c 2000 2>/dev/null # core file size (512-byte blocks) +ulimit -c 0 2>/dev/null # core file size (512-byte blocks) ulimit -v unlimited 2>/dev/null # maximum virtual memory size ulimit -m unlimited 2>/dev/null # maximum physical memory size # disable output in case of scp if [ ! -z "$PS1" ]; then echo - echo Düsseldorf Advanced Materials Simulation Kit - DAMASK - echo Max-Planck-Institut für Eisenforschung, Düsseldorf + echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK + echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf echo https://damask.mpie.de echo echo Using environment with ... @@ -64,14 +61,21 @@ if [ ! -z "$PS1" ]; then echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" if [ "x$PETSC_DIR" != "x" ]; then echo "PETSc location $PETSC_DIR" - [[ `python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"` == $PETSC_DIR ]] \ - || echo " ~~> "`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"` + [[ $(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") == $PETSC_DIR ]] \ + || echo " ~~> "$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") fi - [[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH" + [[ "x$PETSC_ARCH" == "x" ]] \ + || echo "PETSc architecture $PETSC_ARCH" echo "MSC.Marc/Mentat $MSC_ROOT" echo - echo -n "heap size/MiB "; echo "`ulimit -d`/1024" | bc - echo -n "stack size/MiB "; echo "`ulimit -s`/1024" | bc + echo -n "heap size " + [[ "$(ulimit -d)" == "unlimited" ]] \ + && echo "unlimited" \ + || echo $(python -c "import math; size=$(( 1024*$(ulimit -d) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + echo -n "stack size " + [[ "$(ulimit -s)" == "unlimited" ]] \ + && echo "unlimited" \ + || echo $(python -c "import math; size=$(( 1024*$(ulimit -s) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") fi export DAMASK_NUM_THREADS diff --git a/DAMASK_env.zsh b/DAMASK_env.zsh index e434a97cf..082fa8007 100644 --- a/DAMASK_env.zsh +++ b/DAMASK_env.zsh @@ -1,12 +1,7 @@ # sets up an environment for DAMASK on zsh # usage: source DAMASK_env.zsh - -if [ "$OSTYPE" = "linux-gnu" ] || [ "$OSTYPE" = 'linux' ]; then - DAMASK_ROOT=$(readlink -f "`dirname ${(%):-%N}`") -else - print 'Not done yet' -fi +DAMASK_ROOT=${0:a:h} # defining set() allows to source the same file for tcsh and zsh, with and without space around = set() { @@ -15,45 +10,36 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set -# if DAMASK_BIN is present and not in $PATH, add it +# add DAMASK_BIN if present but not yet in $PATH MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:` if [[ ( "x$DAMASK_BIN" != "x" ) && ( "x$MATCH" = "x" ) ]]; then export PATH=$DAMASK_BIN:$PATH fi SOLVER=`which DAMASK_spectral 2>/dev/null` -if [ "x$SOLVER" = "x" ]; then - export SOLVER='Not found!' -fi PROCESSING=`which postResults 2>/dev/null` -if [ "x$PROCESSING" = "x" ]; then - export PROCESSING='Not found!' -fi if [ "x$DAMASK_NUM_THREADS" = "x" ]; then DAMASK_NUM_THREADS=1 fi # according to http://software.intel.com/en-us/forums/topic/501500 # this seems to make sense for the stack size -FREE=`which free 2>/dev/null` -if [ "x$FREE" != "x" ]; then +if [ "`which free 2>/dev/null`" != "free not found" ]; then freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'` - heap=`expr $freeMem / 2` - stack=`expr $freeMem / 2` # http://superuser.com/questions/220059/what-parameters-has-ulimit - ulimit -s $stack 2>/dev/null # maximum stack size (kB) - ulimit -d $heap 2>/dev/null # maximum heap size (kB) + ulimit -s `expr $freeMem / $DAMASK_NUM_THREADS / 2` 2>/dev/null # maximum stack size (kB) + ulimit -d `expr $freeMem / 2` 2>/dev/null # maximum heap size (kB) fi -ulimit -c 2000 2>/dev/null # core file size (512-byte blocks) +ulimit -c 0 2>/dev/null # core file size (512-byte blocks) ulimit -v unlimited 2>/dev/null # maximum virtual memory size ulimit -m unlimited 2>/dev/null # maximum physical memory size # disable output in case of scp if [ ! -z "$PS1" ]; then echo - echo Düsseldorf Advanced Materials Simulation Kit - DAMASK - echo Max-Planck-Institut für Eisenforschung, Düsseldorf + echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK + echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf echo https://damask.mpie.de echo echo Using environment with ... @@ -63,13 +49,21 @@ if [ ! -z "$PS1" ]; then echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS" if [ "x$PETSC_DIR" != "x" ]; then echo "PETSc location $PETSC_DIR" - [[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR` + [[ $(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") == $PETSC_DIR ]] \ + || echo " ~~> "$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") fi - [[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH" + [[ "x$PETSC_ARCH" == "x" ]] \ + || echo "PETSc architecture $PETSC_ARCH" echo "MSC.Marc/Mentat $MSC_ROOT" echo - echo -n "heap size/MiB "; echo "`ulimit -d`/1024" | bc - echo -n "stack size/MiB "; echo "`ulimit -s`/1024" | bc + echo -n "heap size " + [[ "$(ulimit -d)" == "unlimited" ]] \ + && echo "unlimited" \ + || echo $(python -c "import math; size=$(( 1024*$(ulimit -d) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + echo -n "stack size " + [[ "$(ulimit -s)" == "unlimited" ]] \ + && echo "unlimited" \ + || echo $(python -c "import math; size=$(( 1024*$(ulimit -s) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") fi export DAMASK_NUM_THREADS From 5e1cac2e5e6456794149b3a9d62a7e6f036369ed Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 29 Jul 2016 17:46:22 -0400 Subject: [PATCH 11/15] pretty-printing of CONFIG --- CONFIG | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/CONFIG b/CONFIG index b701ec3c7..9be3a95fb 100644 --- a/CONFIG +++ b/CONFIG @@ -1,10 +1,11 @@ # "set"-syntax needed only for tcsh (but works with bash and zsh) # DAMASK_ROOT will be expanded -set DAMASK_BIN=${DAMASK_ROOT}/bin -set DAMASK_NUM_THREADS = 4 +set DAMASK_BIN = ${DAMASK_ROOT}/bin -set MSC_ROOT=/opt/MSC -set MARC_VERSION=2015 +set DAMASK_NUM_THREADS = 4 -set ABAQUS_VERSION=6.14-5 +set MSC_ROOT = /opt/MSC +set MARC_VERSION = 2015 + +set ABAQUS_VERSION = 6.14-5 From afff0b86143c534f480c34bd5c35081eca19629d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 29 Jul 2016 17:48:40 -0400 Subject: [PATCH 12/15] fixed STDOUT error no files on command line translate to "filename" being empty list. Cannot test for filename[0] then... --- processing/post/perceptualUniformColorMap.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/processing/post/perceptualUniformColorMap.py b/processing/post/perceptualUniformColorMap.py index 85acd073b..96e676337 100755 --- a/processing/post/perceptualUniformColorMap.py +++ b/processing/post/perceptualUniformColorMap.py @@ -49,7 +49,7 @@ parser.set_defaults(right = (0.0,0.0,0.0)) (options,filename) = parser.parse_args() if options.format not in outtypes: - parser.error('invalid format: "%s" (can be %s).'%(options.format,', '.join(outtypes))) + parser.error('invalid format: "{}" (choices: {}).'.format(options.format,', '.join(outtypes))) if options.N < 2: parser.error('too few steps (need at least 2).') @@ -59,10 +59,9 @@ if options.trim[0] < -1.0 or \ options.trim[0] >= options.trim[1]: parser.error('invalid trim range (-1 +1).') - -name = options.format if filename[0] is None\ +name = options.format if filename == [] \ else filename[0] -output = sys.stdout if filename[0] is None\ +output = sys.stdout if filename == [] \ else open(os.path.basename(filename[0])+extensions[outtypes.index(options.format)],'w') colorLeft = damask.Color(options.colormodel.upper(), list(options.left)) From 0bbf54e0e46fbb021661b3b15ab4260dcf3cea39 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 29 Jul 2016 17:49:29 -0400 Subject: [PATCH 13/15] switched to string.format() method --- processing/post/addCauchy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addCauchy.py b/processing/post/addCauchy.py index 2f69cb043..2d366b198 100755 --- a/processing/post/addCauchy.py +++ b/processing/post/addCauchy.py @@ -68,7 +68,7 @@ for name in filenames: # ------------------------------------------ assemble header -------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.labels_append(['%i_Cauchy'%(i+1) for i in xrange(9)]) # extend ASCII header with new labels + table.labels_append(['{}_Cauchy'.format(i+1) for i in xrange(9)]) # extend ASCII header with new labels table.head_write() # ------------------------------------------ process data ------------------------------------------ From 1503cc19fdd346276caa2787e34c6c5ff41931d3 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 29 Jul 2016 20:40:40 -0400 Subject: [PATCH 14/15] fixed excessive line length --- DAMASK_env.sh | 12 ++++++++++-- DAMASK_env.zsh | 12 ++++++++++-- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/DAMASK_env.sh b/DAMASK_env.sh index d9b5b004f..da033afbe 100644 --- a/DAMASK_env.sh +++ b/DAMASK_env.sh @@ -71,11 +71,19 @@ if [ ! -z "$PS1" ]; then echo -n "heap size " [[ "$(ulimit -d)" == "unlimited" ]] \ && echo "unlimited" \ - || echo $(python -c "import math; size=$(( 1024*$(ulimit -d) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + || echo $(python -c \ + "import math; \ + size=$(( 1024*$(ulimit -d) )); \ + print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ + ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])") echo -n "stack size " [[ "$(ulimit -s)" == "unlimited" ]] \ && echo "unlimited" \ - || echo $(python -c "import math; size=$(( 1024*$(ulimit -s) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + || echo $(python -c \ + "import math; \ + size=$(( 1024*$(ulimit -s) )); \ + print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ + ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])") fi export DAMASK_NUM_THREADS diff --git a/DAMASK_env.zsh b/DAMASK_env.zsh index 082fa8007..fa9a33987 100644 --- a/DAMASK_env.zsh +++ b/DAMASK_env.zsh @@ -59,11 +59,19 @@ if [ ! -z "$PS1" ]; then echo -n "heap size " [[ "$(ulimit -d)" == "unlimited" ]] \ && echo "unlimited" \ - || echo $(python -c "import math; size=$(( 1024*$(ulimit -d) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + || echo $(python -c \ + "import math; \ + size=$(( 1024*$(ulimit -d) )); \ + print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ + ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])") echo -n "stack size " [[ "$(ulimit -s)" == "unlimited" ]] \ && echo "unlimited" \ - || echo $(python -c "import math; size=$(( 1024*$(ulimit -s) )); print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), ['bytes', 'KiB', 'MiB', 'GiB', 'TiB', 'EiB', 'ZiB'][int(math.log(size,2) / 10) if size else 0])") + || echo $(python -c \ + "import math; \ + size=$(( 1024*$(ulimit -s) )); \ + print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \ + ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])") fi export DAMASK_NUM_THREADS From 304fdf1ebe47da15eba1e85e46ee786eced6d194 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Fri, 29 Jul 2016 20:41:15 -0400 Subject: [PATCH 15/15] can deal with "veterans" and "newbies" meaning over ride existing with new --- processing/post/addCalculation.py | 80 ++++++++++++++----------------- 1 file changed, 37 insertions(+), 43 deletions(-) diff --git a/processing/post/addCalculation.py b/processing/post/addCalculation.py index 5e8c155b1..1d603e9fb 100755 --- a/processing/post/addCalculation.py +++ b/processing/post/addCalculation.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python2 # -*- coding: UTF-8 no BOM -*- import os,re,sys @@ -6,10 +6,15 @@ import math import numpy as np from optparse import OptionParser import damask +import collections scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) +def listify(x): + return (x if isinstance(x, collections.Iterable) else [x]) + + # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -55,13 +60,9 @@ for i in xrange(len(options.formulas)): if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - buffered = False) - output = damask.ASCIItable(name = name, - buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name, + buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------- @@ -129,53 +130,46 @@ for name in filenames: while outputAlive and table.data_read(): # read next data line of ASCII table specials['_row_'] += 1 # count row - output.data_clear() # ------------------------------------------ calculate one result to get length of labels --------- if firstLine: firstLine = False - labelDim = {} - for label in [x for x in options.labels]: - labelDim[label] = np.size(eval(evaluator[label])) - if labelDim[label] == 0: options.labels.remove(label) + resultDim = {} + for label in list(options.labels): # iterate over stable copy + resultDim[label] = np.size(eval(evaluator[label])) # get dimension of formula[label] + if resultDim[label] == 0: options.labels.remove(label) # remove label if invalid result + + veterans = list(set(options.labels)&set(table.labels(raw=False)+table.labels(raw=True)) ) # intersection of requested and existing + newbies = list(set(options.labels)-set(table.labels(raw=False)+table.labels(raw=True)) ) # requested but not existing + + for veteran in list(veterans): + if resultDim[veteran] != table.label_dimension(veteran): + damask.util.croak('{} is ignored due to inconsistent dimension...'.format(veteran)) + veterans.remove(veteran) # ignore culprit + for newby in newbies: + table.labels_append(['{}_{}'.format(i+1,newby) for i in xrange(resultDim[newby])] + if resultDim[newby] > 1 else newby) # ------------------------------------------ assemble header --------------------------------------- - output.labels_clear() - tabLabels = table.labels() - for label in tabLabels: - dim = labelDim[label] if label in options.labels \ - else table.label_dimension(label) - output.labels_append(['{}_{}'.format(i+1,label) for i in xrange(dim)] if dim > 1 else label) - - for label in options.labels: - if label in tabLabels: continue - output.labels_append(['{}_{}'.format(i+1,label) for i in xrange(labelDim[label])] - if labelDim[label] > 1 - else label) - - output.info = table.info - output.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - output.head_write() + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + table.head_write() # ------------------------------------------ process data ------------------------------------------ - for label in output.labels(): - oldIndices = table.label_indexrange(label) - Nold = max(1,len(oldIndices)) # Nold could be zero for new columns - Nnew = len(output.label_indexrange(label)) - output.data_append(eval(evaluator[label]) if label in options.labels and - (options.condition is None or eval(eval(evaluator_condition))) - else np.tile([table.data[i] for i in oldIndices] - if label in tabLabels - else np.nan, - np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns + if options.condition is None or eval(eval(evaluator_condition)): # condition for veteran replacement fulfilled + for veteran in veterans: # evaluate formulae that overwrite + table.data[table.label_index(veteran): + table.label_index(veteran)+table.label_dimension(veteran)] = \ + listify(eval(evaluator[veteran])) + + for newby in newbies: # evaluate formulae that append + table.data_append(listify(eval(evaluator[newby]))) - outputAlive = output.data_write() # output processed line + outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- - table.input_close() # close ASCII tables - output.close() # close ASCII tables - + table.close() # close ASCII table +