diff --git a/PRIVATE b/PRIVATE index 64432754c..9dc7065be 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 +Subproject commit 9dc7065beedf2a097aa60a656bbfa52e55d7147c diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 681691c9e..6f3c91fff 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -25,44 +25,44 @@ module constitutive use kinematics_cleavage_opening use kinematics_slipplane_opening use kinematics_thermal_expansion - + implicit none private - + integer, public, protected :: & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState - + interface module subroutine plastic_none_init end subroutine plastic_none_init - + module subroutine plastic_isotropic_init end subroutine plastic_isotropic_init - + module subroutine plastic_phenopowerlaw_init end subroutine plastic_phenopowerlaw_init - + module subroutine plastic_kinehardening_init end subroutine plastic_kinehardening_init - + module subroutine plastic_dislotwin_init end subroutine plastic_dislotwin_init module subroutine plastic_disloUCLA_init end subroutine plastic_disloUCLA_init - + module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init - - + + 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) :: & @@ -75,20 +75,20 @@ module constitutive 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 end subroutine plastic_phenopowerlaw_LpAndItsTangent - + pure module subroutine plastic_kinehardening_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) :: & @@ -101,7 +101,7 @@ module constitutive 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 real(pReal), intent(in) :: & @@ -110,13 +110,13 @@ module constitutive instance, & of end subroutine plastic_dislotwin_LpAndItsTangent - + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,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 real(pReal), intent(in) :: & @@ -125,14 +125,14 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_LpAndItsTangent - + module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Mp, Temperature, volume, ip, el) 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 real(pReal), intent(in) :: & @@ -149,15 +149,15 @@ module constitutive 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 end subroutine plastic_isotropic_LiAndItsTangent - - + + module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -165,7 +165,7 @@ module constitutive instance, & of end subroutine plastic_isotropic_dotState - + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -173,7 +173,7 @@ module constitutive instance, & of end subroutine plastic_phenopowerlaw_dotState - + module subroutine plastic_kinehardening_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -181,7 +181,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_dotState - + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -201,8 +201,8 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) integer, intent(in) :: & ip, & !< current integration point @@ -213,11 +213,11 @@ module constitutive real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< deformation gradient Fp !< plastic deformation gradient end subroutine plastic_nonlocal_dotState - + module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & instance, & @@ -225,19 +225,19 @@ module constitutive real(pReal), intent(in) :: & T end subroutine plastic_dislotwin_dependentState - + module subroutine plastic_disloUCLA_dependentState(instance,of) integer, intent(in) :: & instance, & of end subroutine plastic_disloUCLA_dependentState - - module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + + module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp end subroutine plastic_nonlocal_dependentState @@ -249,7 +249,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_deltaState - + module subroutine plastic_nonlocal_deltaState(Mp,ip,el) integer, intent(in) :: & ip, & @@ -267,48 +267,48 @@ module constitutive ip, & !< integration point el !< element end function plastic_dislotwin_homogenizedC - - module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + + module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - - + + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_isotropic_results - + module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_phenopowerlaw_results - + module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results - + module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_dislotwin_results - + module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_disloUCLA_results - + module subroutine plastic_nonlocal_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results - + end interface - + public :: & plastic_nonlocal_updateCompatibility, & constitutive_init, & @@ -321,7 +321,7 @@ module constitutive constitutive_collectDotState, & constitutive_collectDeltaState, & constitutive_results - + contains @@ -355,18 +355,18 @@ subroutine constitutive_init if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init - + !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init - + write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) - + constitutive_plasticity_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0 - + PhaseLoop2:do ph = 1,material_Nphase !-------------------------------------------------------------------------------------------------- ! partition and inititalize state @@ -398,7 +398,7 @@ function constitutive_homogenizedC(ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) @@ -412,23 +412,23 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient integer :: & ho, & !< homogenization tme, & !< thermal member position instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -439,7 +439,7 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloUCLA_dependentState(instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (Fe,Fp,ip,el) + call plastic_nonlocal_dependentState (F,Fp,ip,el) end select plasticityType end subroutine constitutive_microstructure @@ -451,7 +451,7 @@ end subroutine constitutive_microstructure ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) + S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -473,49 +473,49 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & tme !< thermal member position integer :: & i, j, instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + Mp = matmul(matmul(transpose(Fi),Fi),S) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - + case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + case (PLASTICITY_ISOTROPIC_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_KINEHARDENING_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) - + case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, & temperature(ho)%p(tme),geometry_plastic_nonlocal_IPvolume0(ip,el),ip,el) - + case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + case (PLASTICITY_DISLOUCLA_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + end select plasticityType - + do i=1,3; do j=1,3 dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) @@ -545,7 +545,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dS, & !< derivative of Li with respect to S dLi_dFi - + real(pReal), dimension(3,3) :: & my_Li, & !< intermediate velocity gradient FiInv, & @@ -557,11 +557,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & integer :: & k, i, j, & instance, of - + Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_isotropic_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -571,10 +571,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal end select plasticityType - + Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el)) kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType @@ -590,12 +590,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop - + FiInv = math_inv33(Fi) detFi = math_det33(Fi) Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = matmul(FiInv,Li) - + do i = 1,3; do j = 1,3 dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) @@ -621,10 +621,10 @@ pure function constitutive_initialFi(ipc, ip, el) integer :: & phase, & homog, offset - + constitutive_initialFi = math_I3 phase = material_phaseAt(ipc,el) - + KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType @@ -640,7 +640,7 @@ end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic/intermediate deformation gradients depending on the selected elastic law +!> the elastic/intermediate deformation gradients depending on the selected elastic law !! (so far no case switch because only Hooke is implemented) !-------------------------------------------------------------------------------------------------- subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) @@ -690,17 +690,17 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & d !< counter in degradation loop integer :: & i, j - + ho = material_homogenizationAt(el) C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) - + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 end select degradationType enddo DegradationLoop - + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration @@ -715,7 +715,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) +subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -724,7 +724,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FeArray, & !< elastic deformation gradient + FArray, & !< elastic deformation gradient FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient @@ -771,7 +771,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme), & subdt,ip,el) end select plasticityType @@ -858,32 +858,32 @@ subroutine constitutive_results do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) call results_closeGroup(results_addGroup(group)) - + group = trim(group)//'/plastic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) select case(phase_plasticity(p)) - + case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) - + call plastic_isotropic_results(phase_plasticityInstance(p),group) + case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) - + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + call plastic_kinehardening_results(phase_plasticityInstance(p),group) case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) - + call plastic_dislotwin_results(phase_plasticityInstance(p),group) + case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group) - + call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) + call plastic_nonlocal_results(phase_plasticityInstance(p),group) end select - - enddo + + enddo end subroutine constitutive_results diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 42b8a34f5..13d629e2b 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -11,7 +11,7 @@ submodule(constitutive) plastic_nonlocal IPvolume => geometry_plastic_nonlocal_IPvolume0, & IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 - + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -44,11 +44,11 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:), allocatable :: & totalNslip !< total number of active slip systems for each instance !END DEPRECATED - + real(pReal), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between me and my neighbors - - enum, bind(c) + + enum, bind(c) enumerator :: & undefined_ID, & rho_sgl_mob_edg_pos_ID, & @@ -73,7 +73,7 @@ submodule(constitutive) plastic_nonlocal v_scr_neg_ID, & gamma_ID end enum - + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume @@ -135,27 +135,27 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID !< ID of each post result output - + end type tParameters - + type :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & - tau_pass, & - tau_Back + tau_pass, & + tau_Back end type tNonlocalMicrostructure - + type :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & - rhoSglMobile, & ! iRhoU + rhoSglMobile, & ! iRhoU rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_neg, & rho_sgl_mob_scr_pos, & @@ -176,14 +176,15 @@ submodule(constitutive) plastic_nonlocal v_scr_pos, & v_scr_neg end type tNonlocalState - + type(tNonlocalState), allocatable, dimension(:) :: & deltaState, & dotState, & - state - + state, & + state0 + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - + type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure contains @@ -210,8 +211,8 @@ module subroutine plastic_nonlocal_init extmsg = '', & structure character(len=pStringLen), dimension(:), allocatable :: outputs - integer :: NofMyPhase - + integer :: NofMyPhase + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' @@ -226,6 +227,7 @@ module subroutine plastic_nonlocal_init allocate(param(maxNinstances)) allocate(state(maxNinstances)) + allocate(state0(maxNinstances)) allocate(dotState(maxNinstances)) allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) @@ -238,13 +240,14 @@ module subroutine plastic_nonlocal_init associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & + st0 => state0(phase_plasticityInstance(p)), & del => deltaState(phase_plasticityInstance(p)), & dst => microstructure(phase_plasticityInstance(p)), & config => config_phase(p)) prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=0.0_pReal) - + structure = config%getString('lattice_structure') ! This data is read in already in lattice @@ -293,20 +296,20 @@ module subroutine plastic_nonlocal_init prm%colinearSystem(s1) = s2 enddo enddo - + prm%rhoSglEdgePos0 = config%getFloats('rhosgledgepos0', requiredSize=size(prm%Nslip)) prm%rhoSglEdgeNeg0 = config%getFloats('rhosgledgeneg0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewPos0 = config%getFloats('rhosglscrewpos0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewNeg0 = config%getFloats('rhosglscrewneg0', requiredSize=size(prm%Nslip)) prm%rhoDipEdge0 = config%getFloats('rhodipedge0', requiredSize=size(prm%Nslip)) prm%rhoDipScrew0 = config%getFloats('rhodipscrew0', requiredSize=size(prm%Nslip)) - + prm%lambda0 = config%getFloats('lambda0', requiredSize=size(prm%Nslip)) prm%burgers = config%getFloats('burgers', requiredSize=size(prm%Nslip)) - + prm%lambda0 = math_expand(prm%lambda0,prm%Nslip) prm%burgers = math_expand(prm%burgers,prm%Nslip) - + prm%minDipoleHeight_edge = config%getFloats('minimumdipoleheightedge', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_screw = config%getFloats('minimumdipoleheightscrew', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge,prm%Nslip) @@ -314,7 +317,7 @@ module subroutine plastic_nonlocal_init allocate(prm%minDipoleHeight(prm%totalNslip,2)) prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw - + prm%peierlsstress_edge = config%getFloats('peierlsstressedge', requiredSize=size(prm%Nslip)) prm%peierlsstress_screw = config%getFloats('peierlsstressscrew', requiredSize=size(prm%Nslip)) prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge,prm%Nslip) @@ -326,7 +329,7 @@ module subroutine plastic_nonlocal_init prm%significantRho = config%getFloat('significantrho') prm%significantN = config%getFloat('significantn', 0.0_pReal) prm%CFLfactor = config%getFloat('cflfactor',defaultVal=2.0_pReal) - + prm%atomicVolume = config%getFloat('atomicvolume') prm%Dsd0 = config%getFloat('selfdiffusionprefactor') !,'dsd0' prm%selfDiffusionEnergy = config%getFloat('selfdiffusionenergy') !,'qsd' @@ -336,7 +339,7 @@ module subroutine plastic_nonlocal_init prm%solidSolutionEnergy = config%getFloat('solidsolutionenergy') prm%solidSolutionSize = config%getFloat('solidsolutionsize') prm%solidSolutionConcentration = config%getFloat('solidsolutionconcentration') - + prm%p = config%getFloat('p') prm%q = config%getFloat('q') prm%viscosity = config%getFloat('viscosity') @@ -345,62 +348,62 @@ module subroutine plastic_nonlocal_init ! ToDo: discuss logic prm%rhoSglScatter = config%getFloat('rhosglscatter') prm%rhoSglRandom = config%getFloat('rhosglrandom',0.0_pReal) - if (config%keyExists('/rhosglrandom/')) & + if (config%keyExists('/rhosglrandom/')) & prm%rhoSglRandomBinning = config%getFloat('rhosglrandombinning',0.0_pReal) !ToDo: useful default? ! if (rhoSglRandom(instance) < 0.0_pReal) & ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) & - + prm%surfaceTransmissivity = config%getFloat('surfacetransmissivity',defaultVal=1.0_pReal) prm%grainboundaryTransmissivity = config%getFloat('grainboundarytransmissivity',defaultVal=-1.0_pReal) prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' lambda0' - + if (any(prm%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgePos0' if (any(prm%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgeNeg0' if (any(prm%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewPos0' if (any(prm%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewNeg0' if (any(prm%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipEdge0' if (any(prm%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipScrew0' - + if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' peierlsstress' if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' minDipoleHeight' - + if (prm%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' viscosity' if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' selfDiffusionEnergy' if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' fattack' if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' doublekinkwidth' if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' Dsd0' if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' atomicVolume' ! ToDo: in disloUCLA/dislotwin, the atomic volume is given as a factor - + if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' significantN' if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' significantrho' if (prm%atolshear <= 0.0_pReal) extmsg = trim(extmsg)//' atolshear' if (prm%atolrho <= 0.0_pReal) extmsg = trim(extmsg)//' atolrho' if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' CFLfactor' - - if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' - if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' - + + if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' + if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) & extmsg = trim(extmsg)//' linetensionEffect' if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) & extmsg = trim(extmsg)//' edgeJogFactor' - + if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionEnergy' if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionSize' if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionConcentration' - if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' + if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & - extmsg = trim(extmsg)//' surfaceTransmissivity' - + extmsg = trim(extmsg)//' surfaceTransmissivity' + if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & - extmsg = trim(extmsg)//' fEdgeMultiplication' + extmsg = trim(extmsg)//' fEdgeMultiplication' endif slipActive @@ -470,39 +473,40 @@ module subroutine plastic_nonlocal_init 'gamma ' ]) * prm%totalNslip !< "basic" microstructural state variables that are independent from other state variables sizeDependentState = size([ 'rhoForest ']) * prm%totalNslip !< microstructural state variables that depend on other state variables sizeState = sizeDotState + sizeDependentState & - + size([ 'velocityEdgePos ','velocityEdgeNeg ', & + + size([ 'velocityEdgePos ','velocityEdgeNeg ', & 'velocityScrewPos ','velocityScrewNeg ', & 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%totalNslip !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - + call material_allocatePlasticState(p,NofMyPhase,sizeState,sizeDotState,sizeDeltaState) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - + totalNslip(phase_plasticityInstance(p)) = prm%totalNslip - + + st0%rho => plasticState(p)%state0 (0*prm%totalNslip+1:10*prm%totalNslip,:) stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:) del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:) plasticState(p)%aTolState(1:10*prm%totalNslip) = prm%aTolRho - + stt%rhoSgl => plasticState(p)%state (0*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSgl => plasticState(p)%dotState (0*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSgl => plasticState(p)%deltaState (0*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoSglMobile => plasticState(p)%state (0*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rhoSglMobile => plasticState(p)%dotState (0*prm%totalNslip+1: 4*prm%totalNslip,:) del%rhoSglMobile => plasticState(p)%deltaState (0*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%totalNslip+1: 1*prm%totalNslip,:) dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%totalNslip+1: 1*prm%totalNslip,:) del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%totalNslip+1: 1*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%totalNslip+1: 2*prm%totalNslip,:) dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%totalNslip+1: 2*prm%totalNslip,:) del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%totalNslip+1: 2*prm%totalNslip,:) - + stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%totalNslip+1: 3*prm%totalNslip,:) dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%totalNslip+1: 3*prm%totalNslip,:) del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%totalNslip+1: 3*prm%totalNslip,:) @@ -510,46 +514,46 @@ module subroutine plastic_nonlocal_init stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%totalNslip+1: 4*prm%totalNslip,:) del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rhoSglImmobile => plasticState(p)%state (4*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%totalNslip+1: 5*prm%totalNslip,:) dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%totalNslip+1: 5*prm%totalNslip,:) del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%totalNslip+1: 5*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%totalNslip+1: 6*prm%totalNslip,:) dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%totalNslip+1: 6*prm%totalNslip,:) del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%totalNslip+1: 6*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%totalNslip+1: 7*prm%totalNslip,:) dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%totalNslip+1: 7*prm%totalNslip,:) del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%totalNslip+1: 7*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%totalNslip+1: 8*prm%totalNslip,:) del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoDip => plasticState(p)%state (8*prm%totalNslip+1:10*prm%totalNslip,:) dot%rhoDip => plasticState(p)%dotState (8*prm%totalNslip+1:10*prm%totalNslip,:) del%rhoDip => plasticState(p)%deltaState (8*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%rho_dip_edg => plasticState(p)%state (8*prm%totalNslip+1: 9*prm%totalNslip,:) dot%rho_dip_edg => plasticState(p)%dotState (8*prm%totalNslip+1: 9*prm%totalNslip,:) del%rho_dip_edg => plasticState(p)%deltaState (8*prm%totalNslip+1: 9*prm%totalNslip,:) - + stt%rho_dip_scr => plasticState(p)%state (9*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho_dip_scr => plasticState(p)%dotState (9*prm%totalNslip+1:10*prm%totalNslip,:) del%rho_dip_scr => plasticState(p)%deltaState (9*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%gamma => plasticState(p)%state (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) dot%gamma => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) del%gamma => plasticState(p)%deltaState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%aTolState(10*prm%totalNslip + 1:11*prm%totalNslip ) = prm%aTolShear plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%accumulatedSlip => plasticState(p)%state(10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) - + stt%rho_forest => plasticState(p)%state (11*prm%totalNslip + 1:12*prm%totalNslip ,1:NofMyPhase) stt%v => plasticState(p)%state (12*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase) stt%v_edg_pos => plasticState(p)%state (12*prm%totalNslip + 1:13*prm%totalNslip ,1:NofMyPhase) @@ -565,10 +569,10 @@ module subroutine plastic_nonlocal_init plasticState(p)%state0 = plasticState(p)%state enddo - + allocate(compatibility(2,maxval(totalNslip),maxval(totalNslip),nIPneighbors,& discretization_nIP,discretization_nElem), source=0.0_pReal) - + ! BEGIN DEPRECATED---------------------------------------------------------------------------------- allocate(iRhoU(maxval(totalNslip),4,maxNinstances), source=0) allocate(iRhoB(maxval(totalNslip),4,maxNinstances), source=0) @@ -601,7 +605,7 @@ module subroutine plastic_nonlocal_init iRhoD(s,c,phase_plasticityInstance(p)) = l enddo enddo - + l = l + param(phase_plasticityInstance(p))%totalNslip ! shear(rates) l = l + param(phase_plasticityInstance(p))%totalNslip ! rho_forest @@ -619,20 +623,20 @@ module subroutine plastic_nonlocal_init enddo if (iD(param(phase_plasticityInstance(p))%totalNslip,2,phase_plasticityInstance(p)) /= plasticState(p)%sizeState) & call IO_error(0, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')') - - + + endif myPhase2 - + enddo initializeInstances ! END DEPRECATED------------------------------------------------------------------------------------ - + contains !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- subroutine stateInit(phase,NofMyPhase) - + integer,intent(in) ::& phase, & NofMyPhase @@ -647,7 +651,7 @@ module subroutine plastic_nonlocal_init phasemember real(pReal), dimension(2) :: & noise, & - rnd + rnd real(pReal) :: & meanDensity, & totalVolume, & @@ -655,14 +659,14 @@ module subroutine plastic_nonlocal_init minimumIpVolume real(pReal), dimension(NofMyPhase) :: & volume - - + + instance = phase_plasticityInstance(phase) associate(prm => param(instance), stt => state(instance)) - - ! randomly distribute dislocation segments on random slip system and of random type in the volume + + ! randomly distribute dislocation segments on random slip system and of random type in the volume if (prm%rhoSglRandom > 0.0_pReal) then - + ! get the total volume of the instance do e = 1,discretization_nElem do i = 1,discretization_nIP @@ -672,7 +676,7 @@ module subroutine plastic_nonlocal_init totalVolume = sum(volume) minimumIPVolume = minval(volume) densityBinning = prm%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) - + ! subsequently fill random ips with dislocation segments until we reach the desired overall density meanDensity = 0.0_pReal do while(meanDensity < prm%rhoSglRandom) @@ -701,9 +705,9 @@ module subroutine plastic_nonlocal_init enddo enddo endif - + end associate - + end subroutine stateInit end subroutine plastic_nonlocal_init @@ -712,15 +716,15 @@ end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) - +module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) + integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp - + integer :: & ph, & !< phase of, & !< offset @@ -761,11 +765,12 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & - rho_neighbor + rho0, & + rho_neighbor0 real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), & totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor @@ -774,30 +779,30 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) neighbor_rhoTotal ! total density at neighboring material point real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & m ! direction of dislocation motion - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) - + ns = prm%totalNslip - + rho = getRho(instance,of,ip,el) - + stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & - + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) - - - ! coefficients are corrected for the line tension effect + + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) + + + ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) - if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc + if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then do s = 1,ns correction = ( 1.0_pReal - prm%linetensionEffect & + prm%linetensionEffect & * log(0.35_pReal * prm%burgers(s) * sqrt(max(stt%rho_forest(s,of),prm%significantRho))) & / log(0.35_pReal * prm%burgers(s) * 1e6_pReal)) ** 2.0_pReal - myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) + myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) enddo else myInteractionMatrix = prm%interactionSlipSlip @@ -815,21 +820,21 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - + rho0 = getRho0(instance,of,ip,el) if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then - invFe = math_inv33(Fe) invFp = math_inv33(Fp) - - rho_edg_delta = rho(:,mob_edg_pos) - rho(:,mob_edg_neg) - rho_scr_delta = rho(:,mob_scr_pos) - rho(:,mob_scr_neg) - + invFe = matmul(Fp,math_inv33(F)) + + rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) + rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) + rhoExcess(1,1:ns) = rho_edg_delta rhoExcess(2,1:ns) = rho_scr_delta - + FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) - + !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - + nRealNeighbors = 0.0_pReal neighbor_rhoTotal = 0.0_pReal do n = 1,nIPneighbors @@ -839,16 +844,16 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) if (neighbor_instance == instance) then - + nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor = getRho(instance,no,neighbor_ip,neighbor_el) - - rho_edg_delta_neighbor(:,n) = rho_neighbor(:,mob_edg_pos) - rho_neighbor(:,mob_edg_neg) - rho_scr_delta_neighbor(:,n) = rho_neighbor(:,mob_scr_pos) - rho_neighbor(:,mob_scr_neg) - - neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor(:,edg)),2) - neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) - + rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el) + + rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) + rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) + + neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) + neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) + connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - discretization_IPcoords(1:3,el+neighbor_ip-1)) normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) @@ -867,18 +872,18 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta_neighbor(:,n) = rho_scr_delta endif enddo - + neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor - + !* loop through the slip systems and calculate the dislocation gradient by !* 1. interpolation of the excess density in the neighorhood !* 2. interpolation of the dead dislocation density in the central volume m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_transverse - + do s = 1,ns - + ! gradient from interpolation of neighboring excess density ... do c = 1,2 do dir = 1,3 @@ -891,27 +896,27 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) enddo invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') - + rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences)) enddo - + ! ... plus gradient from deads ... rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize - + ! ... normalized with the total density ... rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pReal + nRealNeighbors) rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pReal + nRealNeighbors) - + rhoExcessGradient_over_rho = 0.0_pReal where(rhoTotal > 0.0_pReal) & rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal - + ! ... gives the local stress correction when multiplied with a factor dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * pi) & * (rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) - + enddo endif @@ -932,7 +937,7 @@ end subroutine plastic_nonlocal_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief calculates kinetics +!> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tauThreshold, c, Temperature, instance, of) @@ -945,17 +950,17 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress - + real(pReal), dimension(param(instance)%totalNslip), intent(out) :: & v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) - + integer :: & ns, & !< short notation for the total number of active slip systems s !< index of my current slip system - real(pReal) :: & - tauRel_P, & + real(pReal) :: & + tauRel_P, & tauRel_S, & tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers @@ -976,22 +981,22 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & criticalStress_P, & !< maximum obstacle strength criticalStress_S, & !< maximum obstacle strength mobility !< dislocation mobility - + associate(prm => param(instance)) ns = prm%totalNslip v = 0.0_pReal dv_dtau = 0.0_pReal dv_dtauNS = 0.0_pReal - - + + if (Temperature > 0.0_pReal) then do s = 1,ns if (abs(tau(s)) > tauThreshold(s)) then - + !* Peierls contribution !* Effective stress includes non Schmid constributions !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive meanfreepath_P = prm%burgers(s) jumpWidth_P = prm%burgers(s) @@ -1006,15 +1011,15 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & if (tauEff < criticalStress_P) then dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) & - * tauRel_P**(prm%p-1.0_pReal) + * tauRel_P**(prm%p-1.0_pReal) else dtPeierls_dtau = 0.0_pReal endif - - + + !* Contribution from solid solution strengthening !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = abs(tau(s)) - tauThreshold(s) meanfreepath_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) jumpWidth_S = prm%solidSolutionSize * prm%burgers(s) @@ -1030,32 +1035,32 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q & * activationVolume_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal) & - * tauRel_S**(prm%p-1.0_pReal) + * tauRel_S**(prm%p-1.0_pReal) else dtSolidSolution_dtau = 0.0_pReal endif - - + + !* viscous glide velocity - + tauEff = abs(tau(s)) - tauThreshold(s) mobility = prm%burgers(s) / prm%viscosity vViscous = mobility * tauEff - - - !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of - !* free flight at glide velocity in between. + + + !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of + !* free flight at glide velocity in between. !* adopt sign from resolved stress - + v(s) = sign(1.0_pReal,tau(s)) & / (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) dv_dtau(s) = v(s) * v(s) * (dtSolidSolution_dtau / meanfreepath_S & + mobility / (vViscous * vViscous)) - dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P + dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P endif enddo endif - + #ifdef DEBUGTODO write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal @@ -1089,8 +1094,8 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to Tstar (9x9 matrix) - - + + integer :: & instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems @@ -1114,20 +1119,20 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - + !*** shortcut for mapping ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) - + instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance)) ns = prm%totalNslip - - !*** shortcut to state variables + + !*** shortcut to state variables rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) - - + + !*** get resolved shear stress !*** for screws possible non-schmid contributions are also taken into account do s = 1,ns @@ -1145,18 +1150,18 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & forall (t = 1:4) & tauNS(1:ns,t) = tauNS(1:ns,t) + dst%tau_back(:,of) tau = tau + dst%tau_back(:,of) - - + + !*** get dislocation velocity and its tangent and store the velocity in the state array - - ! edges + + ! edges call plastic_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & tau(1:ns), tauNS(1:ns,1), dst%tau_pass(1:ns,of), & 1, Temperature, instance, of) v(1:ns,2) = v(1:ns,1) dv_dtau(1:ns,2) = dv_dtau(1:ns,1) dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) - + !screws if (size(prm%nonSchmidCoeff) == 0) then forall(t = 3:4) @@ -1173,17 +1178,17 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & endif stt%v(:,of) = pack(v,.true.) - + !*** Bauschinger effect forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) - - + + gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * prm%burgers(1:ns) - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + do s = 1,ns Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1191,13 +1196,13 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) & + prm%Schmid(i,j,s) & - * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) enddo - - + + end associate - + end subroutine plastic_nonlocal_LpAndItsTangent @@ -1205,7 +1210,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_deltaState(Mp,ip,el) - + integer, intent(in) :: & ip, & el @@ -1234,23 +1239,23 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance)) ns = totalNslip(instance) - - !*** shortcut to state variables + + !*** shortcut to state variables forall (s = 1:ns, t = 1:4) & v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) forall (s = 1:ns, c = 1:2) & dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of) - + rho = getRho(instance,of,ip,el) rhoDip = rho(:,dip) - + !**************************************************************************** !*** dislocation remobilization (bauschinger effect) where(rho(:,imm) * v < 0.0_pReal) @@ -1261,48 +1266,48 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) elsewhere deltaRhoRemobilization(:,mob) = 0.0_pReal deltaRhoRemobilization(:,imm) = 0.0_pReal - endwhere + endwhere deltaRhoRemobilization(:,dip) = 0.0_pReal - + !**************************************************************************** !*** calculate dipole formation and dissociation by stress change - + !*** calculate limits for stable dipole height do s = 1,prm%totalNslip tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - + dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) - + dUpper = max(dUpper,prm%minDipoleHeight) deltaDUpper = dUpper - dUpperOld - - + + !*** dissociation by stress increase deltaRhoDipole2SingleStress = 0.0_pReal forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) - + forall (t=1:4) & deltaRhoDipole2SingleStress(1:ns,t) = -0.5_pReal & * deltaRhoDipole2SingleStress(1:ns,(t-1)/2+9) - + forall (s = 1:ns, c = 1:2) & - plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) - + plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) + plasticState(ph)%deltaState(:,of) = 0.0_pReal del%rho(:,of) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1320,9 +1325,9 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & +module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) - + integer, intent(in) :: & ip, & !< current integration point el !< current element number @@ -1332,11 +1337,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient - + integer :: & - ph, & + ph, & instance, & !< current instance of this plasticity neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems @@ -1358,6 +1363,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & + rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotFlux, & !< density evolution by flux @@ -1366,17 +1372,17 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) - my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & v, & !< current dislocation glide velocity - my_v, & !< dislocation glide velocity of central ip - neighbor_v, & !< dislocation glide velocity of enighboring ip + v0, & + neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws @@ -1399,19 +1405,19 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion - + logical :: & considerEnteringFlux, & considerLeavingFlux - + p = material_phaseAt(1,el) o = material_phasememberAt(1,ip,el) - + if (timestep <= 0.0_pReal) then plasticState(p)%dotState = 0.0_pReal return endif - + ph = material_phaseAt(1,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance), & @@ -1419,25 +1425,27 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & dot => dotState(instance), & stt => state(instance)) ns = totalNslip(instance) - + tau = 0.0_pReal gdot = 0.0_pReal - - rho = getRho(instance,o,ip,el) + + rho = getRho(instance,o,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - + rho0 = getRho0(instance,o,ip,el) + my_rhoSgl0 = rho0(:,sgl) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(p)%state(iV (s,t,instance),o) endforall - - + + !**************************************************************************** !*** Calculate shear rate - + forall (t = 1:4) & gdot(1:ns,t) = rhoSgl(1:ns,t) * prm%burgers(1:ns) * v(1:ns,t) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1446,29 +1454,29 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot endif #endif - - - + + + !**************************************************************************** !*** calculate limits for stable dipole height - + do s = 1,ns ! loop over slip systems tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,o) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - + dLower = prm%minDipoleHeight dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) dUpper = max(dUpper,dLower) - + !**************************************************************************** !*** calculate dislocation multiplication rhoDotMultiplication = 0.0_pReal @@ -1481,29 +1489,32 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & * sqrt(stt%rho_forest(s,o)) / prm%lambda0(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall - + else isBCC rhoDotMultiplication(1:ns,1:4) = spread( & (sum(abs(gdot(1:ns,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(1:ns,3:4)),2)) & * sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4) endif isBCC - - + + forall (s = 1:ns, t = 1:4) + v0(s,t) = plasticState(p)%state0(iV(s,t,instance),o) + endforall + !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) rhoDotFlux = 0.0_pReal if (.not. phase_localPlasticity(material_phaseAt(1,el))) then - + !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%CFLfactor * abs(v) * timestep & + .and. prm%CFLfactor * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & - maxval(abs(v), abs(gdot) > 0.0_pReal & - .and. prm%CFLfactor * abs(v) * timestep & + maxval(abs(v0), abs(gdot) > 0.0_pReal & + .and. prm%CFLfactor * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep write(6,'(a)') '<< CONST >> enforcing cutback !!!' @@ -1512,70 +1523,69 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback return endif - - + + !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! !*** opposite sign to our p vector in the (s,p,n) triplet !!! - + m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_direction m(1:3,1:ns,3) = -prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse - - my_Fe = Fe(1:3,1:3,1,ip,el) - my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el)) - + + my_F = F(1:3,1:3,1,ip,el) + my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) + neighbors: do n = 1,nIPneighbors - + neighbor_el = IPneighborhood(1,n,ip,el) neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_n = IPneighborhood(3,n,ip,el) np = material_phaseAt(1,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el) - + opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) - + if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) - neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el)) + neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) + neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F endif - - + + !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity - !* (also nonlocal constitutive law with local properties) that is at least a little bit + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) that is at least a little bit !* compatible. !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. - !* The entering flux from my neighbor will be distributed on my slip systems according to the + !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility - + considerEnteringFlux = .false. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below - neighbor_rhoSgl = 0.0_pReal + neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below if (neighbor_n > 0) then if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & considerEnteringFlux = .true. endif - + enteringFlux: if (considerEnteringFlux) then forall (s = 1:ns, t = 1:4) - neighbor_v(s,t) = plasticState(np)%state(iV (s,t,neighbor_instance),no) - neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), & + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no), & 0.0_pReal) endforall - - where (neighbor_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & - .or. neighbor_rhoSgl < prm%significantRho) & - neighbor_rhoSgl = 0.0_pReal + + where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & + .or. neighbor_rhoSgl0 < prm%significantRho) & + neighbor_rhoSgl0 = 0.0_pReal normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & @@ -1586,9 +1596,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & do t = 1,4 c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me - .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & + if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me + .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & @@ -1602,107 +1612,104 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & enddo enddo endif enteringFlux - - + + !* FLUX FROM ME TO MY NEIGHBOR - !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). + !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. !* So the net flux in the direction of my neighbor is equal to zero: !* leaving flux to neighbor == entering flux from opposite neighbor !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. - + considerLeavingFlux = .true. if (opposite_n > 0) then if (phase_plasticity(material_phaseAt(1,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & considerLeavingFlux = .false. endif - + leavingFlux: if (considerLeavingFlux) then - my_rhoSgl = rhoSgl - my_v = v - normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)), & + * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration area = IParea(n,ip,el) * norm2(normal_me2neighbor) - normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length + normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 c = (t + 1) / 2 - if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) - if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density + if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal endif - lineLength = my_rhoSgl(s,t) * my_v(s,t) & + lineLength = my_rhoSgl0(s,t) * v0(s,t) & * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & - * sign(1.0_pReal, my_v(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point endif enddo enddo endif leavingFlux - + enddo neighbors endif - - - + + + !**************************************************************************** !*** calculate dipole formation and annihilation - + !*** formation by glide - + do c = 1,2 rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns)& * rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) & - rhoDotSingle2DipoleGlide(1:ns,2*c) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4)) enddo - - + + !*** athermal annihilation - + rhoDotAthermalAnnihilation = 0.0_pReal - - forall (c=1:2) & + + forall (c=1:2) & rhoDotAthermalAnnihilation(1:ns,c+8) = -2.0_pReal * dLower(1:ns,c) / prm%burgers(1:ns) & * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system - + if (lattice_structure(ph) == LATTICE_fcc_ID) & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & * 0.25_pReal * sqrt(stt%rho_forest(s,o)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor - - - + + + !*** thermally activated annihilation of edge dipoles by climb - + rhoDotThermalAnnihilation = 0.0_pReal selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (KB * Temperature)) vClimb = prm%atomicVolume * selfDiffusion / ( KB * Temperature ) & @@ -1713,12 +1720,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = 0.0_pReal rhoDot = rhoDotFlux & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation + + rhoDotThermalAnnihilation #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & @@ -1747,7 +1753,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & if ( any(rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%aTolRho) & .or. any(rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%aTolRho)) then #ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif @@ -1765,21 +1771,21 @@ end subroutine plastic_nonlocal_dotState !-------------------------------------------------------------------------------------------------- !> @brief Compatibility update -!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip +!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) - +module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation ! crystal orientation - + integer :: & Nneighbors, & ! number of neighbors - n, & ! neighbor index + n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor ph, & @@ -1792,13 +1798,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) s2 ! slip system index (my neighbor) real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& - nIPneighbors) :: & - my_compatibility ! my_compatibility for current element and ip + nIPneighbors) :: & + my_compatibility ! my_compatibility for current element and ip real(pReal) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & + logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & belowThreshold type(rotation) :: mis @@ -1808,32 +1814,32 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) instance = phase_plasticityInstance(ph) ns = totalNslip(instance) associate(prm => param(instance)) - + !*** start out fully compatible my_compatibility = 0.0_pReal - - forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal - + + forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal + !*** Loop thrugh neighbors and check whether there is any compatibility. - + neighbors: do n = 1,Nneighbors neighbor_e = IPneighborhood(1,n,i,e) neighbor_i = IPneighborhood(2,n,i,e) - - + + !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config - + if (neighbor_e <= 0 .or. neighbor_i <= 0) then forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = sqrt(prm%surfaceTransmissivity) cycle endif - - + + !* PHASE BOUNDARY - !* If we encounter a different nonlocal phase at the neighbor, + !* If we encounter a different nonlocal phase at the neighbor, !* we consider this to be a real "physical" phase boundary, so completely incompatible. - !* If one of the two phases has a local plasticity law, + !* If one of the two phases has a local plasticity law, !* we do not consider this to be a phase boundary, so completely compatible. neighbor_phase = material_phaseAt(1,neighbor_e) if (neighbor_phase /= ph) then @@ -1841,8 +1847,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal cycle endif - - + + !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) if (prm%grainboundaryTransmissivity >= 0.0_pReal) then @@ -1854,16 +1860,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) endif cycle endif - - + + !* GRAIN BOUNDARY ? !* Compatibility defined by relative orientation of slip systems: !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) @@ -1878,7 +1884,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) * abs(math_inner(prm%slip_direction(1:3,s1), & mis%rotate(prm%slip_direction(1:3,s2)))) enddo neighborSlipSystems - + my_compatibilitySum = 0.0_pReal belowThreshold = .true. do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns))) @@ -1896,33 +1902,33 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal enddo mySlipSystems endif - + enddo neighbors - + compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility - + end associate - + end subroutine plastic_nonlocal_updateCompatibility - + !-------------------------------------------------------------------------------------------------- !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- function getRho(instance,of,ip,el) - + integer, intent(in) :: instance, of,ip,el real(pReal), dimension(param(instance)%totalNslip,10) :: getRho - + associate(prm => param(instance)) getRho = reshape(state(instance)%rho(:,of),[prm%totalNslip,10]) - + ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - + where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & getRho = 0.0_pReal @@ -1931,6 +1937,31 @@ function getRho(instance,of,ip,el) end function getRho +!-------------------------------------------------------------------------------------------------- +!> @brief returns copy of current dislocation densities from state +!> @details raw values is rectified +!-------------------------------------------------------------------------------------------------- +function getRho0(instance,of,ip,el) + + integer, intent(in) :: instance, of,ip,el + real(pReal), dimension(param(instance)%totalNslip,10) :: getRho0 + + associate(prm => param(instance)) + + getRho0 = reshape(state0(instance)%rho(:,of),[prm%totalNslip,10]) + + ! ensure positive densities (not for imm, they have a sign) + getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) + getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) + + where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & + getRho0 = 0.0_pReal + + end associate + +end function getRho0 + + !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8356d0a65..83ef3fefc 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -22,10 +22,10 @@ module crystallite use discretization use lattice use results - - implicit none + + implicit none private - + real(pReal), dimension(:,:,:), allocatable, public :: & crystallite_dt !< requested time increment of each grain real(pReal), dimension(:,:,:), allocatable :: & @@ -33,7 +33,7 @@ module crystallite crystallite_subFrac, & !< already calculated fraction of increment crystallite_subStep !< size of next integration step type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation + crystallite_orientation !< current orientation real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_P !< 1st Piola-Kirchhoff stress per grain @@ -74,33 +74,33 @@ module crystallite crystallite_converged, & !< convergence flag crystallite_todo, & !< flag to indicate need for further computation crystallite_localPlasticity !< indicates this grain to have purely local constitutive law - + type :: tOutput !< new requested output (per phase) character(len=pStringLen), allocatable, dimension(:) :: & label end type tOutput type(tOutput), allocatable, dimension(:) :: output_constituent - + type :: tNumerics integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress !< stress loop limit + nStress !< stress loop limit real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback subStepSizeLp, & !< size of first substep when cutback in Lp calculation subStepSizeLi, & !< size of first substep when cutback in Li calculation stepIncreaseCryst, & !< increase of next substep size when previous substep converged - rTol_crystalliteState, & !< relative tolerance in state loop + rTol_crystalliteState, & !< relative tolerance in state loop rTol_crystalliteStress, & !< relative tolerance in stress loop aTol_crystalliteStress !< absolute tolerance in stress loop end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? - + procedure(), pointer :: integrateState - + public :: & crystallite_init, & crystallite_stress, & @@ -116,7 +116,7 @@ contains !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- subroutine crystallite_init - + logical, dimension(discretization_nIP,discretization_nElem) :: devNull integer :: & c, & !< counter in integration point component loop @@ -126,13 +126,13 @@ subroutine crystallite_init iMax, & !< maximum number of integration points eMax, & !< maximum number of elements myNcomponents !< number of components at current IP - + write(6,'(/,a)') ' <<<+- crystallite init -+>>>' - + cMax = homogenization_maxNgrains iMax = discretization_nIP eMax = discretization_nElem - + allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -172,23 +172,23 @@ subroutine crystallite_init allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - + num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal) num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal) - + num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal) num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal) num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal) num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal) num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal) - + num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) - + num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nStress = config_numerics%getInt ('nstress', defaultVal=40) - + if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') @@ -199,12 +199,12 @@ subroutine crystallite_init if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState') if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress') if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress') - + if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - + if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - + select case(numerics_integrator) case(1) integrateState => integrateStateFPI @@ -217,11 +217,11 @@ subroutine crystallite_init case(5) integrateState => integrateStateRKCK45 end select - + allocate(output_constituent(size(config_phase))) do c = 1, size(config_phase) #if defined(__GFORTRAN__) - allocate(output_constituent(c)%label(1)) + allocate(output_constituent(c)%label(1)) output_constituent(c)%label(1)= 'GfortranBug86277' output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label ) if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::] @@ -250,28 +250,28 @@ subroutine crystallite_init enddo; enddo enddo !$OMP END PARALLEL DO - + if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? - + crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedF0 = crystallite_F0 crystallite_partionedF = crystallite_F0 - + call crystallite_orientations() - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) - call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fp(1:3,1:3,c,i,e), & + call constitutive_microstructure(crystallite_partionedF0(1:3,1:3,c,i,e), & + crystallite_partionedFp0(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo enddo enddo !$OMP END PARALLEL DO - + devNull = crystallite_stress() call crystallite_stressTangent @@ -283,7 +283,7 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif - + call debug_info call debug_reset #endif @@ -295,7 +295,7 @@ end subroutine crystallite_init !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - + logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress real(pReal), intent(in), optional :: & dummyArgumentToPreventInternalCompilerErrorWithGCC @@ -308,7 +308,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s - + #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & .and. FEsolving_execElem(1) <= debug_e & @@ -494,7 +494,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_stress = .false. elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) + crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) enddo enddo elementLooping5 @@ -675,12 +675,12 @@ end subroutine crystallite_stressTangent !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- subroutine crystallite_orientations - + integer & c, & !< counter in integration point component loop i, & !< counter in integration point loop e !< counter in element loop - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -688,7 +688,7 @@ subroutine crystallite_orientations call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) enddo; enddo; enddo !$OMP END PARALLEL DO - + nonlocalPresent: if (any(plasticState%nonLocal)) then !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -706,7 +706,7 @@ end subroutine crystallite_orientations !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(ipc,ip,el, tensor33) - + real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3) :: T @@ -714,7 +714,7 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33) el, & ip, & ipc - + T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) @@ -731,11 +731,11 @@ subroutine crystallite_results real(pReal), allocatable, dimension(:,:,:) :: selected_tensors type(rotation), allocatable, dimension(:) :: selected_rotations character(len=256) :: group,lattice_label - + do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) do o = 1, size(output_constituent(p)%label) select case (output_constituent(p)%label(o)) @@ -792,19 +792,19 @@ subroutine crystallite_results end select enddo enddo - + contains !------------------------------------------------------------------------------------------------ !> @brief select tensors for output !------------------------------------------------------------------------------------------------ function select_tensors(dataset,instance) - + integer, intent(in) :: instance real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset real(pReal), allocatable, dimension(:,:,:) :: select_tensors integer :: e,i,c,j - + allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -818,20 +818,20 @@ subroutine crystallite_results enddo enddo enddo - + end function select_tensors - - + + !-------------------------------------------------------------------------------------------------- !> @brief select rotations for output -!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- function select_rotations(dataset,instance) - + integer, intent(in) :: instance type(rotation), dimension(:,:,:), intent(in) :: dataset type(rotation), allocatable, dimension(:) :: select_rotations integer :: e,i,c,j - + allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -845,7 +845,7 @@ subroutine crystallite_results enddo enddo enddo - + end function select_rotations end subroutine crystallite_results @@ -856,12 +856,12 @@ end subroutine crystallite_results !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- logical function integrateStress(ipc,ip,el,timeFraction) - + integer, intent(in):: el, & ! element index ip, & ! integration point index ipc ! grain index real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - + real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep @@ -916,7 +916,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) jacoCounterLi ! counters to check for Jacobian update external :: & dgesv - + !* be pessimistic integrateStress = .false. #ifdef DEBUG @@ -938,7 +938,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess Liguess_old = Liguess - + invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el)) failedInversionFp: if (all(dEq0(invFp_current))) then #ifdef DEBUG @@ -964,13 +964,13 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif failedInversionFi - + !* start Li loop with normal step length NiterationStressLi = 0 jacoCounterLi = 0 steplengthLi = 1.0_pReal residuumLi_old = 0.0_pReal - + LiLoop: do NiterationStressLi = NiterationStressLi + 1 LiLoopLimit: if (NiterationStressLi > num%nStress) then @@ -981,18 +981,18 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif LiLoopLimit - + invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) Fi_new = math_inv33(invFi_new) detInvFi = math_det33(invFi_new) - + !* start Lp loop with normal step length NiterationStressLp = 0 jacoCounterLp = 0 steplengthLp = 1.0_pReal residuumLp_old = 0.0_pReal Lpguess_old = Lpguess - + LpLoop: do NiterationStressLp = NiterationStressLp + 1 LpLoopLimit: if (NiterationStressLp > num%nStress) then @@ -1003,18 +1003,18 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif LpLoopLimit - + !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law - + B = math_I3 - dt*Lpguess Fe = matmul(matmul(A,B), invFi_new) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration - + !* calculate plastic velocity gradient and its tangent from constitutive law call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) - + #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -1032,7 +1032,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error num%aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive - + if (any(IEEE_is_NaN(residuumLp))) then #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & @@ -1160,7 +1160,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif cycle LiLoop endif - + !* calculate Jacobian for correction term if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then temp_33 = matmul(matmul(A,B),invFi_current) @@ -1197,11 +1197,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif - + deltaLi = - math_9to33(work) endif jacoCounterLi = jacoCounterLi + 1 - + Liguess = Liguess + steplengthLi * deltaLi #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & @@ -1211,7 +1211,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) endif #endif enddo LiLoop - + !* calculate new plastic and elastic deformation gradient invFp_new = matmul(invFp_current,B) invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize @@ -1302,7 +1302,7 @@ subroutine integrateStateFPI #endif ! store previousDotState and previousDotState2 - + !$OMP PARALLEL DO PRIVATE(p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1329,7 +1329,7 @@ subroutine integrateStateFPI call update_dependentState call update_stress(1.0_pReal) call update_dotState(1.0_pReal) - + !$OMP PARALLEL !$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1342,7 +1342,7 @@ subroutine integrateStateFPI zeta = damper(plasticState(p)%dotState (:,c), & plasticState(p)%previousDotState (:,c), & plasticState(p)%previousDotState2(:,c)) - + residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - ( plasticState(p)%dotState (:,c) * zeta & @@ -1350,18 +1350,18 @@ subroutine integrateStateFPI ) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - residuum_plastic(1:sizeDotState) + - residuum_plastic(1:sizeDotState) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + zeta = damper(sourceState(p)%p(s)%dotState (:,c), & sourceState(p)%p(s)%previousDotState (:,c), & sourceState(p)%p(s)%previousDotState2(:,c)) @@ -1432,12 +1432,12 @@ subroutine integrateStateFPI !> @brief calculate the damping for correction of state and dot state !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(current,previous,previous2) - + real(pReal), dimension(:), intent(in) ::& current, previous, previous2 - + real(pReal) :: dot_prod12, dot_prod22 - + dot_prod12 = dot_product(current - previous, previous - previous2) dot_prod22 = dot_product(previous - previous2, previous - previous2) if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then @@ -1445,7 +1445,7 @@ subroutine integrateStateFPI else damper = 1.0_pReal endif - + end function damper end subroutine integrateStateFPI @@ -1480,7 +1480,7 @@ subroutine integrateStateAdaptiveEuler c, & s, & sizeDotState - + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & @@ -1501,14 +1501,14 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) plasticState(p)%state(1:sizeDotState,c) = & plasticState(p)%state(1:sizeDotState,c) + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) sourceState(p)%p(s)%state(1:sizeDotState,c) = & @@ -1530,17 +1530,17 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = & residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) @@ -1549,13 +1549,13 @@ subroutine integrateStateAdaptiveEuler sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo - + endif enddo; enddo; enddo !$OMP END PARALLEL DO if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - + end subroutine integrateStateAdaptiveEuler @@ -1720,23 +1720,23 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) - + residuum_plastic(1:sizeDotState,g,i,e) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & ! why transpose? Better to transpose constant DB * crystallite_subdt(g,i,e) - + plasticState(p)%dotState(:,cc) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) - + residuum_source(1:sizeDotState,s,g,i,e) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) @@ -1744,13 +1744,13 @@ subroutine integrateStateRKCK45 sourceState(p)%p(s)%dotState(:,cc) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B) enddo - + endif enddo; enddo; enddo !$OMP END PARALLEL DO call update_state(1.0_pReal) - + ! --- relative residui and state convergence --- !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) @@ -1759,16 +1759,16 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,cc), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + crystallite_todo(g,i,e) = & crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,cc), & @@ -1783,7 +1783,7 @@ subroutine integrateStateRKCK45 call update_stress(1.0_pReal) call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - + end subroutine integrateStateRKCK45 @@ -1792,7 +1792,7 @@ end subroutine integrateStateRKCK45 !> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - + if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... where( .not. crystallite_localPlasticity) crystallite_converged = .false. @@ -1810,7 +1810,7 @@ subroutine setConvergenceFlag e, & !< element index in element loop i, & !< integration point index in ip loop g !< grain index in grain loop - + !OMP DO PARALLEL PRIVATE do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1826,7 +1826,7 @@ end subroutine setConvergenceFlag !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,aTol) - + real(pReal), intent(in), dimension(:) ::& residuum, state, aTol real(pReal) :: & @@ -1949,11 +1949,11 @@ subroutine update_dotState(timeFraction) g, & !< grain index in grain loop p, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE (p,c,NaN) @@ -1963,9 +1963,9 @@ subroutine update_dotState(timeFraction) !$OMP FLUSH(nonlocalStop) if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe, & + crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_Fp, & + crystallite_partionedFp0, & crystallite_subdt(g,i,e)*timeFraction, g,i,e) p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) @@ -1980,7 +1980,7 @@ subroutine update_dotState(timeFraction) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity + if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity end subroutine update_DotState @@ -1995,11 +1995,11 @@ subroutine update_deltaState mySize, & myOffset, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN) @@ -2016,23 +2016,23 @@ subroutine update_deltaState myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c))) - + if (.not. NaN) then - + plasticState(p)%state(myOffset + 1: myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1: myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) do s = 1, phase_Nsources(p) myOffset = sourceState(p)%p(s)%offsetDeltaState mySize = sourceState(p)%p(s)%sizeDeltaState NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%deltaState(1:mySize,c))) - + if (.not. NaN) then sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) = & sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) + sourceState(p)%p(s)%deltaState(1:mySize,c) endif enddo endif - + crystallite_todo(g,i,e) = .not. NaN if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken crystallite_converged(g,i,e) = .false. @@ -2042,7 +2042,7 @@ subroutine update_deltaState enddo; enddo; enddo !$OMP END PARALLEL DO if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity - + end subroutine update_deltaState