From f97800658f768e5a26b9235e31a0cc58ff434fd0 Mon Sep 17 00:00:00 2001 From: Eureka Pai Date: Thu, 24 May 2018 17:56:09 -0400 Subject: [PATCH 1/6] replaced param(instance) with p => pointer --- src/plastic_isotropic.f90 | 147 +++++++++++++++++++++----------------- 1 file changed, 81 insertions(+), 66 deletions(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index a3a1d5caf..631e9661a 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -53,7 +53,7 @@ module plastic_isotropic dilatation = .false. end type - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable, target, private :: param !< containers of constitutive parameters (len Ninstance) type, private :: tIsotropicState !< internal state aliases real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance @@ -130,6 +130,7 @@ subroutine plastic_isotropic_init(fileUnit) implicit none integer(pInt), intent(in) :: fileUnit + type tParameters, pointer :: p integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: & @@ -185,13 +186,14 @@ subroutine plastic_isotropic_init(fileUnit) if (IO_getTag(line,'[',']') /= '') then ! next section phase = phase + 1_pInt ! advance section counter if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then - instance = phase_plasticityInstance(phase) ! count instances of my constitutive law - allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output + p => param(phase_plasticityInstance(phase)) ! shorthand pointer to parameter object of my constitutive law + allocate(p%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output endif cycle ! skip to next line endif if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase + p => param(instance) chunkPos = IO_stringPos(line) tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key @@ -201,58 +203,58 @@ subroutine plastic_isotropic_init(fileUnit) select case(outputtag) case ('flowstress') plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - param(instance)%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID + p%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag case ('strainrate') plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - param(instance)%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID + p%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag end select case ('/dilatation/') - param(instance)%dilatation = .true. + p%dilatation = .true. case ('tau0') - param(instance)%tau0 = IO_floatValue(line,chunkPos,2_pInt) + p%tau0 = IO_floatValue(line,chunkPos,2_pInt) case ('gdot0') - param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) + p%gdot0 = IO_floatValue(line,chunkPos,2_pInt) case ('n') - param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) + p%n = IO_floatValue(line,chunkPos,2_pInt) case ('h0') - param(instance)%h0 = IO_floatValue(line,chunkPos,2_pInt) + p%h0 = IO_floatValue(line,chunkPos,2_pInt) case ('h0_slope','slopelnrate') - param(instance)%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt) + p%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt) case ('tausat') - param(instance)%tausat = IO_floatValue(line,chunkPos,2_pInt) + p%tausat = IO_floatValue(line,chunkPos,2_pInt) case ('tausat_sinhfita') - param(instance)%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt) + p%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt) case ('tausat_sinhfitb') - param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt) + p%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt) case ('tausat_sinhfitc') - param(instance)%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt) + p%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt) case ('tausat_sinhfitd') - param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt) + p%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt) case ('a', 'w0') - param(instance)%a = IO_floatValue(line,chunkPos,2_pInt) + p%a = IO_floatValue(line,chunkPos,2_pInt) case ('taylorfactor') - param(instance)%fTaylor = IO_floatValue(line,chunkPos,2_pInt) + p%fTaylor = IO_floatValue(line,chunkPos,2_pInt) case ('atol_flowstress') - param(instance)%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt) + p%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt) case ('atol_shear') - param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) + p%aTolShear = IO_floatValue(line,chunkPos,2_pInt) case default @@ -269,17 +271,18 @@ subroutine plastic_isotropic_init(fileUnit) myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) instance = phase_plasticityInstance(phase) + p => param(instance) extmsg = '' !-------------------------------------------------------------------------------------------------- ! sanity checks - if (param(instance)%aTolShear <= 0.0_pReal) param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 - if (param(instance)%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0' - if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' - if (param(instance)%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' - if (param(instance)%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat' - if (param(instance)%a <= 0.0_pReal) extmsg = trim(extmsg)//' a' - if (param(instance)%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor' - if (param(instance)%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress' + if (p%aTolShear <= 0.0_pReal) p%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 + if (p%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0' + if (p%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' + if (p%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' + if (p%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat' + if (p%a <= 0.0_pReal) extmsg = trim(extmsg)//' a' + if (p%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor' + if (p%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress' if (extmsg /= '') then extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier call IO_error(211_pInt,ip=instance,ext_msg=extmsg) @@ -287,7 +290,7 @@ subroutine plastic_isotropic_init(fileUnit) !-------------------------------------------------------------------------------------------------- ! Determine size of postResults array outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) - select case(param(instance)%outputID(o)) + select case(p%outputID(o)) case(flowstress_ID,strainrate_ID) mySize = 1_pInt case default @@ -349,13 +352,13 @@ subroutine plastic_isotropic_init(fileUnit) !-------------------------------------------------------------------------------------------------- ! init state - state0(instance)%flowstress = param(instance)%tau0 + state0(instance)%flowstress = p%tau0 state0(instance)%accumulatedShear = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! init absolute state tolerances - stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress - stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear + stateAbsTol(instance)%flowstress = p%aTolFlowstress + stateAbsTol(instance)%accumulatedShear = p%aTolShear endif myPhase enddo initializeInstances @@ -400,6 +403,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) ip, & !< integration point el !< element + type tParameters, pointer :: p + real(pReal), dimension(3,3) :: & Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor real(pReal), dimension(3,3,3,3) :: & @@ -414,7 +419,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - + p => param(instance) + Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) @@ -423,11 +429,11 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) Lp = 0.0_pReal dLp_dTstar99 = 0.0_pReal else - gamma_dot = param(instance)%gdot0 & - * ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) & - **param(instance)%n + gamma_dot = p%gdot0 & + * ( sqrt(1.5_pReal) * norm_Tstar_dev / p%fTaylor / state(instance)%flowstress(of) ) & + **p%n - Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor + Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/p%fTaylor if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -441,13 +447,13 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! Calculation of the tangent of Lp forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & + dLp_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * & Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) & dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal - dLp_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * & + dLp_dTstar99 = math_Plain3333to99(gamma_dot / p%fTaylor * & dLp_dTstar_3333 / norm_Tstar_dev) end if end subroutine plastic_isotropic_LpAndItsTangent @@ -479,9 +485,11 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e ip, & !< integration point el !< element + type tParameters, pointer :: p + real(pReal), dimension(3,3) :: & Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor -real(pReal) :: & + real(pReal) :: & gamma_dot, & !< strainrate norm_Tstar_sph, & !< euclidean norm of Tstar_sph squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph @@ -491,27 +499,28 @@ real(pReal) :: & of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - + p => param(instance) + Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) - if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero - gamma_dot = param(instance)%gdot0 & - * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) & - **param(instance)%n + if (p%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero + gamma_dot = p%gdot0 & + * (sqrt(1.5_pReal) * norm_Tstar_sph / p%fTaylor / state(instance)%flowstress(of) ) & + **p%n - Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor + Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/p%fTaylor !-------------------------------------------------------------------------------------------------- ! Calculation of the tangent of Li forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & + dLi_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * & Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal - dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * & + dLi_dTstar_3333 = gamma_dot / p%fTaylor * & dLi_dTstar_3333 / norm_Tstar_sph else Li = 0.0_pReal @@ -541,6 +550,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element + type tParameters, pointer :: p real(pReal), dimension(6) :: & Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal) :: & @@ -554,10 +564,11 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - + p => param(instance) + !-------------------------------------------------------------------------------------------------- ! norm of (deviatoric) 2nd Piola-Kirchhoff stress - if (param(instance)%dilatation) then + if (p%dilatation) then norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) else Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal @@ -566,26 +577,26 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) end if !-------------------------------------------------------------------------------------------------- ! strain rate - gamma_dot = param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & + gamma_dot = p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & / &!----------------------------------------------------------------------------------- - (param(instance)%fTaylor*state(instance)%flowstress(of) ))**param(instance)%n + (p%fTaylor*state(instance)%flowstress(of) ))**p%n !-------------------------------------------------------------------------------------------------- ! hardening coefficient if (abs(gamma_dot) > 1e-12_pReal) then - if (dEq0(param(instance)%tausat_SinhFitA)) then - saturation = param(instance)%tausat + if (dEq0(p%tausat_SinhFitA)) then + saturation = p%tausat else - saturation = param(instance)%tausat & - + asinh( (gamma_dot / param(instance)%tausat_SinhFitA& - )**(1.0_pReal / param(instance)%tausat_SinhFitD)& - )**(1.0_pReal / param(instance)%tausat_SinhFitC) & - / ( param(instance)%tausat_SinhFitB & - * (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) & + saturation = p%tausat & + + asinh( (gamma_dot / p%tausat_SinhFitA& + )**(1.0_pReal / p%tausat_SinhFitD)& + )**(1.0_pReal / p%tausat_SinhFitC) & + / ( p%tausat_SinhFitB & + * (gamma_dot / p%gdot0)**(1.0_pReal / p%n) & ) endif - hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) & - * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a & + hardening = ( p%h0 + p%h0_slopeLnRate * log(gamma_dot) ) & + * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**p%a & * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation) else hardening = 0.0_pReal @@ -614,6 +625,9 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element + + type tParameters, pointer :: p + real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & plastic_isotropic_postResults @@ -629,10 +643,11 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + p => param(instance) !-------------------------------------------------------------------------------------------------- ! norm of (deviatoric) 2nd Piola-Kirchhoff stress - if (param(instance)%dilatation) then + if (p%dilatation) then norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) else Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal @@ -644,15 +659,15 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) plastic_isotropic_postResults = 0.0_pReal outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) - select case(param(instance)%outputID(o)) + select case(p%outputID(o)) case (flowstress_ID) plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of) c = c + 1_pInt case (strainrate_ID) plastic_isotropic_postResults(c+1_pInt) = & - param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & + p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & / &!---------------------------------------------------------------------------------- - (param(instance)%fTaylor * state(instance)%flowstress(of)) ) ** param(instance)%n + (p%fTaylor * state(instance)%flowstress(of)) ) ** p%n c = c + 1_pInt end select enddo outputsLoop From 982c0fb90ae27665f7d6c403c9b2751b894bf66a Mon Sep 17 00:00:00 2001 From: Eureka Pai Date: Thu, 24 May 2018 18:24:58 -0400 Subject: [PATCH 2/6] replaced param(instance) with p => pointer, corrected error --- src/plastic_isotropic.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 631e9661a..fe415f479 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -130,7 +130,7 @@ subroutine plastic_isotropic_init(fileUnit) implicit none integer(pInt), intent(in) :: fileUnit - type tParameters, pointer :: p + type(tParameters), pointer :: p integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: & From c09a7fd157af44441a899b79401a687d8a36d298 Mon Sep 17 00:00:00 2001 From: Eureka Pai Date: Thu, 24 May 2018 18:31:32 -0400 Subject: [PATCH 3/6] replaced param(instance) with p => pointer, corrected errors --- src/plastic_isotropic.f90 | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index fe415f479..da6d80e60 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -403,7 +403,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) ip, & !< integration point el !< element - type tParameters, pointer :: p + type(tParameters), pointer :: p real(pReal), dimension(3,3) :: & Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor @@ -485,7 +485,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e ip, & !< integration point el !< element - type tParameters, pointer :: p + type(tParameters), pointer :: p real(pReal), dimension(3,3) :: & Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor @@ -526,8 +526,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e Li = 0.0_pReal dLi_dTstar_3333 = 0.0_pReal endif - -end subroutine plastic_isotropic_LiAndItsTangent + end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -550,7 +549,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type tParameters, pointer :: p + type(tParameters), pointer :: p real(pReal), dimension(6) :: & Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal) :: & @@ -626,7 +625,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) ip, & !< integration point el !< element - type tParameters, pointer :: p + type(tParameters), pointer :: p real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & plastic_isotropic_postResults From 8184d51a99ad69ee595889fa50ca8f543a9708c0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 May 2018 09:31:36 +0200 Subject: [PATCH 4/6] following style of more complex constitutive laws offset for different states needs to be computed, so it makes sense to define global and local aliases together. no need to introduce variables for state0 and aTolstate, they are only used once --- src/plastic_isotropic.f90 | 39 ++++++++------------------------------- 1 file changed, 8 insertions(+), 31 deletions(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index da6d80e60..f26eb1ca4 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -61,20 +61,10 @@ module plastic_isotropic accumulatedShear end type - type, private :: tIsotropicAbsTol !< internal alias for abs tolerance in state - real(pReal), pointer :: & ! scalars - flowstress, & - accumulatedShear - end type - type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance state, & - state0, & dotState - type(tIsotropicAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance - stateAbsTol - public :: & plastic_isotropic_init, & plastic_isotropic_LpAndItsTangent, & @@ -263,9 +253,7 @@ subroutine plastic_isotropic_init(fileUnit) enddo parsingFile allocate(state(maxNinstance)) ! internal state aliases - allocate(state0(maxNinstance)) allocate(dotState(maxNinstance)) - allocate(stateAbsTol(maxNinstance)) initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description @@ -334,31 +322,20 @@ subroutine plastic_isotropic_init(fileUnit) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) !-------------------------------------------------------------------------------------------------- -! globally required state aliases - plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase) - plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase) +! locally defined state aliases and initialization of state0 and aTolState -!-------------------------------------------------------------------------------------------------- -! locally defined state aliases state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase) - state0(instance)%flowstress => plasticState(phase)%state0 (1,1:NipcMyPhase) dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase) - stateAbsTol(instance)%flowstress => plasticState(phase)%aTolState(1) + plasticState(phase)%state0(1,1:NipcMyPhase) = p%tau0 + plasticState(phase)%aTolState(1) = p%aTolFlowstress state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase) - state0(instance)%accumulatedShear => plasticState(phase)%state0 (2,1:NipcMyPhase) dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase) - stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) - -!-------------------------------------------------------------------------------------------------- -! init state - state0(instance)%flowstress = p%tau0 - state0(instance)%accumulatedShear = 0.0_pReal - -!-------------------------------------------------------------------------------------------------- -! init absolute state tolerances - stateAbsTol(instance)%flowstress = p%aTolFlowstress - stateAbsTol(instance)%accumulatedShear = p%aTolShear + plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal + plasticState(phase)%aTolState(2) = p%aTolShear + ! global alias + plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase) + plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase) endif myPhase enddo initializeInstances From c7c39922f071f4a31e495b1ba92921d4ad389816 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 May 2018 09:44:14 +0200 Subject: [PATCH 5/6] pointer assignment was done twice pointer is re-assigned automatically, but I found it confusing. Also using automatic left hand side reallocation to simplify handling of outputID --- src/plastic_isotropic.f90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index f26eb1ca4..b30dc441e 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -175,15 +175,11 @@ subroutine plastic_isotropic_init(fileUnit) endif if (IO_getTag(line,'[',']') /= '') then ! next section phase = phase + 1_pInt ! advance section counter - if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then - p => param(phase_plasticityInstance(phase)) ! shorthand pointer to parameter object of my constitutive law - allocate(p%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output - endif cycle ! skip to next line endif if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase - p => param(instance) + p => param(instance) ! shorthand pointer to parameter object of my constitutive law chunkPos = IO_stringPos(line) tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key @@ -193,12 +189,12 @@ subroutine plastic_isotropic_init(fileUnit) select case(outputtag) case ('flowstress') plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - p%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag + p%outputID = [p%outputID,flowstress_ID] case ('strainrate') plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - p%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag + p%outputID = [p%outputID,strainrate_ID] end select case ('/dilatation/') From 2fbe60b949ff7b7fc902a398c9255e16cfee565a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 May 2018 09:54:50 +0200 Subject: [PATCH 6/6] anticipate (proper) change in 23_BasticPETSc_2_PETSc as PRIVATE repository is ahead --- PRIVATE | 2 +- src/spectral_mech_Basic.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PRIVATE b/PRIVATE index 8a2e3ccce..120a0251d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8a2e3ccce34be2a19cdcaecb8456a5dc6c2564f0 +Subproject commit 120a0251d368ee68b35f48e34a249b7eb8bf54a7 diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 8fd2b7d96..7b6931848 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -22,7 +22,7 @@ module spectral_mech_basic private character (len=*), parameter, public :: & - DAMASK_spectral_SolverBasicPETSC_label = 'basicpetsc' + DAMASK_spectral_SolverBasicPETSC_label = 'petsc' !-------------------------------------------------------------------------------------------------- ! derived types