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/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index a3a1d5caf..b30dc441e 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 @@ -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, & @@ -130,6 +120,7 @@ subroutine plastic_isotropic_init(fileUnit) implicit none integer(pInt), intent(in) :: fileUnit + type(tParameters), pointer :: p integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: & @@ -184,14 +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 - 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 - 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) ! shorthand pointer to parameter object of my constitutive law chunkPos = IO_stringPos(line) tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key @@ -201,58 +189,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 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 - param(instance)%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/') - 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 @@ -261,25 +249,24 @@ 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 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 +274,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 @@ -331,31 +318,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 = param(instance)%tau0 - state0(instance)%accumulatedShear = 0.0_pReal - -!-------------------------------------------------------------------------------------------------- -! init absolute state tolerances - stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress - stateAbsTol(instance)%accumulatedShear = param(instance)%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 @@ -400,6 +376,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 +392,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 +402,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 +420,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 +458,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,34 +472,34 @@ 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 dLi_dTstar_3333 = 0.0_pReal endif - -end subroutine plastic_isotropic_LiAndItsTangent + end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -541,6 +522,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 +536,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 +549,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 +597,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 +615,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 +631,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 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