diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index ef8a823d4..beed2112a 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -8,14 +8,7 @@ !! untextured polycrystal !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_isotropic - - enum, bind(c) - enumerator :: & - undefined_ID, & - xi_ID, & - dot_gamma_ID - end enum - + type :: tParameters real(pReal) :: & M, & !< Taylor factor @@ -34,12 +27,12 @@ submodule(constitutive) plastic_isotropic aTol_gamma integer :: & of_debug = 0 - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID logical :: & dilatation + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tIsotropicState real(pReal), pointer, dimension(:) :: & xi, & @@ -56,50 +49,45 @@ submodule(constitutive) plastic_isotropic contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_init integer :: & Ninstance, & - p, i, & + p, & NipcMyPhase, & sizeState, sizeDotState - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' - - write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' - write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' + write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' + Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) & prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) #endif - + prm%xi_0 = config%getFloat('tau0') prm%xi_inf = config%getFloat('tausat') prm%dot_gamma_0 = config%getFloat('gdot0') @@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init prm%a = config%getFloat('a') prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal) prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + prm%dilatation = config%keyExists('/dilatation/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks extmsg = '' @@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m' if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('flowstress') - outputID = xi_ID - case ('strainrate') - outputID = dot_gamma_ID - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['xi ','accumulated_shear']) sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState stt%xi => plasticState(p)%state (1,:) stt%xi = prm%xi_0 dot%xi => plasticState(p)%dotState(1,:) plasticState(p)%aTolState(1) = prm%aTol_xi - + stt%gamma => plasticState(p)%state (2,:) dot%gamma => plasticState(p)%dotState(2,:) plasticState(p)%aTolState(2) = prm%aTol_gamma ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal), dimension(3,3) :: & Mp_dev !< deviatoric part of the Mandel stress real(pReal) :: & @@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev) norm_Mp_dev = sqrt(squarenorm_Mp_dev) - + if (norm_Mp_dev > 0.0_pReal) then dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n - + Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & @@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = 0.0_pReal dLp_dMp = 0.0_pReal end if - + end associate end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & tr !< trace of spherical part of Mandel stress (= 3 x pressure) integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + tr=math_trace33(math_spherical33(Mi)) if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then @@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) Li = 0.0_pReal dLi_dMi = 0.0_pReal endif - + end associate end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_dotState(Mp,instance,of) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + if (prm%dilatation) then norm_Mp = sqrt(math_mul33xx33(Mp,Mp)) else norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp))) endif - + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n - + if (dot_gamma > 1e-12_pReal) then if (dEq0(prm%c_1)) then xi_inf_star = prm%xi_inf @@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) else dot%xi(of) = 0.0_pReal endif - + dot%gamma(of) = dot_gamma ! ToDo: not really used - + end associate end subroutine plastic_isotropic_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_results(instance,group) - integer, intent(in) :: instance + integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (xi_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('flowstress') ! ToDo: should be 'xi' call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa') end select enddo outputsLoop diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 11ce7c494..36cf20af7 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -7,26 +7,13 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_kinehardening - enum, bind(c) - enumerator :: & - undefined_ID, & - crss_ID, & !< critical resolved stress - crss_back_ID, & !< critical resolved back stress - sense_ID, & !< sense of acting shear stress (-1 or +1) - chi0_ID, & !< backstress at last switch of stress sense (positive?) - gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) - accshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - type :: tParameters real(pReal) :: & gdot0, & !< reference shear strain rate for slip n, & !< stress exponent for slip aTolResistance, & aTolShear - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & crss0, & !< initial critical shear stress for slip theta0, & !< initial hardening rate of forward stress for each slip theta1, & !< asymptotic hardening rate of forward stress for each slip @@ -35,21 +22,21 @@ submodule(constitutive) plastic_kinehardening tau1, & tau1_b, & nonSchmidCoeff - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_slipslip !< slip resistance from slip activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & totalNslip, & !< total number of active slip system of_debug = 0 - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance crss, & !< critical resolved stress @@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening gamma0, & !< accumulated shear at last switch of stress sense accshear !< accumulated (absolute) shear end type tKinehardeningState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -72,37 +59,32 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_init integer :: & Ninstance, & - p, i, o, & + p, o, & NipcMyPhase, & sizeState, sizeDeltaState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(deltaState(Ninstance)) - + do p = 1, size(phase_plasticityInstance) if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init dlt => deltaState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)),& config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) then - prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) + prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) endif #endif - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init slipActive: if (prm%totalNslip > 0) then prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) @@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) @@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) - - prm%gdot0 = config%getFloat('gdot0') - prm%n = config%getFloat('n_slip') - + + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n_slip') + ! expand: family => system prm%crss0 = math_expand(prm%crss0, prm%Nslip) prm%tau1 = math_expand(prm%tau1, prm%Nslip) @@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init prm%theta1 = math_expand(prm%theta1, prm%Nslip) prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) - - - + + + !-------------------------------------------------------------------------------------------------- ! sanity checks if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' @@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' - + !ToDo: Any sensible checks for theta? - + endif slipActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance') - outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) - case ('accumulatedshear') - outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) - case ('shearrate') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) - case ('backstress') - outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) - case ('sense') - outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) - case ('chi0') - outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) - case ('gamma0') - outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip sizeState = sizeDotState + sizeDeltaState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init stt%crss = spread(prm%crss0, 2, NipcMyPhase) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%accshear => plasticState(p)%state (startIndex:endIndex,:) @@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + o = plasticState(p)%offsetDeltaState startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%sense => plasticState(p)%state (startIndex :endIndex ,:) dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) @@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - + do i = 1, prm%totalNslip Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo - + end associate end subroutine plastic_kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_dotState(Mp,instance,of) @@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & sumGamma real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg - - + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) dot%accshear(:,of) = abs(gdot_pos+gdot_neg) sumGamma = sum(stt%accshear(:,of)) - - + + dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & * ( prm%theta1 & + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & * exp(-sumGamma*prm%theta0/prm%tau1) & ) - + dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & ( prm%theta1_b + & (prm%theta0_b - prm%theta1_b & + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & ) - + end associate end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates (instantaneous) incremental change of microstructure +!> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & sense - + associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) sense = merge(state(instance)%sense(:,of), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug & @@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) write(6,*) sense,state(instance)%sense(:,of) endif #endif - + !-------------------------------------------------------------------------------------------------- ! switch in sense of shear? where(dNeq(sense,stt%sense(:,of),0.1_pReal)) @@ -402,45 +355,46 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) dlt%chi0 (:,of) = 0.0_pReal dlt%gamma0(:,of) = 0.0_pReal end where - + end associate end subroutine plastic_kinehardening_deltaState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (crss_ID) - call results_writeDataset(group,stt%crss,'xi_sl', & - 'resistance against plastic slip','Pa') + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('resistance') + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', & + 'resistance against plastic slip','Pa') - case(crss_back_ID) - call results_writeDataset(group,stt%crss_back,'tau_back', & - 'back stress against plastic slip','Pa') - - case (sense_ID) - call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') + case('backstress') ! ToDo: should be 'tau_back' + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', & + 'back stress against plastic slip','Pa') + + case ('sense') + if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') + + case ('chi0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') + + case ('gamma0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') + + case ('accumulatedshear') + if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', & + 'plastic shear','1') - case (chi0_ID) - call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') - - case (gamma0_ID) - call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') - - case (accshear_ID) - call results_writeDataset(group,stt%accshear,'gamma_sl', & - 'plastic shear','1') - end select enddo outputsLoop end associate @@ -449,10 +403,11 @@ end subroutine plastic_kinehardening_results !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved +! stress. +!> @details: Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,instance,of, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) @@ -462,44 +417,44 @@ pure subroutine kinetics(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_pos, & gdot_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_pos, & dgdot_dtau_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_pos, & tau_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_pos)) gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where - + where(dNeq0(tau_neg)) gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal end where - + if (present(dgdot_dtau_pos)) then where(dNeq0(gdot_pos)) dgdot_dtau_pos = gdot_pos*prm%n/tau_pos diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 165ad0081..ade8a177f 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -6,19 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_phenopowerlaw - enum, bind(c) - enumerator :: & - undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID - end enum - type :: tParameters real(pReal) :: & gdot0_slip, & !< reference shear strain rate for slip @@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw aTolResistance, & !< absolute tolerance for integration of xi aTolShear, & !< absolute tolerance for integration of gamma aTolTwinfrac !< absolute tolerance for integration of f - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & xi_slip_0, & !< initial critical shear stress for slip xi_twin_0, & !< initial critical shear stress for twin xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & H_int, & !< per family hardening activity (optional) gamma_twin_char !< characteristic shear for twins - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_SlipSlip, & !< slip resistance from slip activity interaction_SlipTwin, & !< slip resistance from twin activity interaction_TwinSlip, & !< twin resistance from slip activity interaction_TwinTwin !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid_slip, & Schmid_twin, & nonSchmid_pos, & @@ -57,13 +44,13 @@ submodule(constitutive) plastic_phenopowerlaw integer :: & totalNslip, & !< total number of active slip system totalNtwin !< total number of active twin systems - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip, & !< number of active slip systems for each family Ntwin !< number of active twin systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_slip, & @@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw gamma_slip, & gamma_twin end type tPhenopowerlawState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -83,7 +70,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_init @@ -91,51 +78,46 @@ module subroutine plastic_phenopowerlaw_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) - + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init slipActive: if (prm%totalNslip > 0) then prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) @@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) - + prm%gdot0_slip = config%getFloat('gdot0_slip') prm%n_slip = config%getFloat('n_slip') prm%a_slip = config%getFloat('a_slip') prm%h0_SlipSlip = config%getFloat('h0_slipslip') - + ! expand: family => system prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) prm%H_int = math_expand(prm%H_int, prm%Nslip) - + ! sanity checks if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' @@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init config%getString('lattice_structure')) prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& config%getFloat('c/a')) - + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - + prm%gdot0_twin = config%getFloat('gdot0_twin') prm%n_twin = config%getFloat('n_twin') prm%spr = config%getFloat('s_pr') prm%h0_TwinTwin = config%getFloat('h0_twintwin') - + ! expand: family => system prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) - + ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' @@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%xi_twin_0(0)) allocate(prm%gamma_twin_char(0)) endif twinActive - + !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then @@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & + size(['tau_twin ','gamma_twin']) * prm%totalNtwin sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNtwin stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) @@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume ! equally (Taylor assumption). Twinning happens only in untwinned volume !-------------------------------------------------------------------------------------------------- @@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & @@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dgdot_dtauslip_pos,dgdot_dtauslip_neg real(pReal), dimension(param(instance)%totalNtwin) :: & gdot_twin,dgdot_dtautwin - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%totalNslip Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) @@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%totalNtwin Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) @@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) enddo twinSystems - + end associate end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) @@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& @@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(param(instance)%totalNslip) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + sumGamma = sum(stt%gamma_slip(:,of)) sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) @@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 - + !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%H_int xi_slip_sat_offset = prm%spr*sqrt(sumF) right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) - + !-------------------------------------------------------------------------------------------------- ! shear rates call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) - + !-------------------------------------------------------------------------------------------------- ! hardening dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) - + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) end associate @@ -437,33 +381,33 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + + case('resistance_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', & + 'resistance against plastic slip','Pa') + case('accumulatedshear_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & + 'plastic shear','1') + + case('resistance_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', & + 'resistance against twinning','Pa') + case('accumulatedshear_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & + 'twinning shear','1') - case (resistance_slip_ID) - call results_writeDataset(group,stt%xi_slip, 'xi_sl', & - 'resistance against plastic slip','Pa') - case (accumulatedshear_slip_ID) - call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & - 'plastic shear','1') - - case (resistance_twin_ID) - call results_writeDataset(group,stt%xi_twin, 'xi_tw', & - 'resistance against twinning','Pa') - case (accumulatedshear_twin_ID) - call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & - 'twinning shear','1') - end select enddo outputsLoop end associate @@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved. +! stress. !> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(Mp,instance,of, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) @@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_slip_pos, & gdot_slip_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_slip_pos, & tau_slip_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_slip_pos)) gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where - + where(dNeq0(tau_slip_neg)) gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where - + if (present(dgdot_dtau_slip_pos)) then where(dNeq0(gdot_slip_pos)) dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos @@ -543,9 +488,9 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. -! twinning is assumed to take place only in untwinned volume. -!> @details Derivates are calculated only optionally. +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. Twinning is assumed to take place only in untwinned volume. +!> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- @@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,& integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & gdot_twin real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & dgdot_dtau_twin - + real(pReal), dimension(param(instance)%totalNtwin) :: & tau_twin integer :: i - + associate(prm => param(instance), stt => state(instance)) - + do i = 1, prm%totalNtwin tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo - + where(tau_twin > 0.0_pReal) gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin else where gdot_twin = 0.0_pReal end where - + if (present(dgdot_dtau_twin)) then where(dNeq0(gdot_twin)) dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin @@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,& dgdot_dtau_twin = 0.0_pReal end where endif - + end associate end subroutine kinetics_twin