From f46d212e47bdc3e675092dc236ee99a0dd3e3b66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 14 Feb 2021 17:29:23 +0100 Subject: [PATCH] simplified access pattern --- src/phase_mechanical.f90 | 23 +++--- src/phase_mechanical_plastic.f90 | 34 ++++----- src/phase_mechanical_plastic_nonlocal.f90 | 85 ++++++++++------------- 3 files changed, 56 insertions(+), 86 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 824c1e257..381820dcc 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -85,11 +85,8 @@ submodule(phase) mechanical logical :: broken end function plastic_dotState - module function plastic_deltaState(co, ip, el, ph, me) result(broken) + module function plastic_deltaState(ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & me logical :: & @@ -114,11 +111,9 @@ submodule(phase) mechanical module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -538,7 +533,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) Fe, Fi_new, co, ip, el) call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, co, ip, el) + S, Fi_new, ph,me) !* update current residuum and check for convergence of loop atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -702,7 +697,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & - r(1:sizeDotState) if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) exit iteration endif @@ -765,7 +760,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -807,7 +802,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -956,7 +951,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D if(broken) return - broken = plastic_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -1329,7 +1324,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & - phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) + phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index cf9faf6c9..c8e9b6271 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -104,7 +104,7 @@ submodule(phase:mechanical) plastic end subroutine dislotungsten_LpAndItsTangent module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me,ip,el) + Mp,Temperature,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -116,9 +116,7 @@ submodule(phase:mechanical) plastic Temperature integer, intent(in) :: & ph, & - me, & - ip, & !< current integration point - el !< current element number + me end subroutine nonlocal_LpAndItsTangent @@ -209,14 +207,12 @@ submodule(phase:mechanical) plastic me end subroutine plastic_kinehardening_deltaState - module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) + module subroutine plastic_nonlocal_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & ph, & - me, & - ip, & - el + me end subroutine plastic_nonlocal_deltaState end interface @@ -244,11 +240,9 @@ end subroutine plastic_init ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -263,14 +257,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & real(pReal), dimension(3,3) :: & Mp !< Mandel stress work conjugate with Lp integer :: & - i, j, me, ph + i, j Mp = matmul(matmul(transpose(Fi),Fi),S) - me = material_phasememberAt(co,ip,el) - ph = material_phaseAt(co,el) - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_NONE_ID) plasticType Lp = 0.0_pReal @@ -286,7 +279,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) case (PLASTICITY_DISLOTWIN_ID) plasticType call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) @@ -391,12 +384,9 @@ end subroutine plastic_dependentState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function plastic_deltaState(co, ip, el, ph, me) result(broken) +module function plastic_deltaState(ph, me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & me logical :: & @@ -421,7 +411,7 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case (PLASTICITY_NONLOCAL_ID) plasticType - call plastic_nonlocal_deltaState(Mp,ph,me,ip,el) + call plastic_nonlocal_deltaState(Mp,ph,me) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case default diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 1f6fce029..b579a3164 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -512,7 +512,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,p,Nconstituents) + if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i) plasticState(p)%state0 = plasticState(p)%state !-------------------------------------------------------------------------------------------------- @@ -576,7 +576,6 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) el integer :: & - ph, & no, & !< neighbor offset neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point @@ -626,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph))) - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -654,7 +653,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me)) @@ -665,7 +664,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) rhoExcess(1,:) = rho_edg_delta rhoExcess(2,:) = rho_scr_delta - FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) + FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal) !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities @@ -680,7 +679,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) if (neighbor_instance == ins(ph)) then nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor0 = getRho0(ph,no,neighbor_ip,neighbor_el) + rho_neighbor0 = getRho0(ph,no) 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) @@ -772,16 +771,14 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me,ip,el) + Mp,Temperature,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp integer, intent(in) :: & ph, & - me, & - ip, & !< current integration point - el !< current element number + me real(pReal), intent(in) :: & Temperature !< temperature @@ -814,7 +811,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) do s = 1,ns @@ -880,18 +877,15 @@ end subroutine nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress integer, intent(in) :: & ph, & - me, & !< offset - ip, & - el + me integer :: & - ph, & !< phase ns, & ! short notation for the total number of active slip systems c, & ! character of dislocation t, & ! type of dislocation @@ -918,7 +912,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me) - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoDip = rho(:,dip) !**************************************************************************** @@ -968,15 +962,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) plasticState(ph)%deltaState(:,me) = 0.0_pReal del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) -#ifdef DEBUG - if (debugConstitutive%extensive & - .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& - .or. .not. debugConstitutive%selective)) then - print'(a,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8) - print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress - endif -#endif - end associate end subroutine plastic_nonlocal_deltaState @@ -1043,10 +1028,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) @@ -1194,7 +1179,6 @@ function rhoDotFlux(timestep,ph,me,ip,el) el !< current element number integer :: & - ph, & neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation @@ -1251,9 +1235,9 @@ function rhoDotFlux(timestep,ph,me,ip,el) gdot = 0.0_pReal - rho = getRho(ph,me,ip,el) + rho = getRho(ph,me) rhoSgl = rho(:,sgl) - rho0 = getRho0(ph,me,ip,el) + rho0 = getRho0(ph,me) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) !ToDo: MD: I think we should use state0 here @@ -1432,7 +1416,6 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor - ph, & neighbor_phase, & ns, & ! number of active slip systems s1, & ! slip system index (me) @@ -1604,7 +1587,7 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nconstituents) +subroutine stateInit(ini,phase,Nconstituents,i) type(tInitialParameters) :: & ini @@ -1631,7 +1614,7 @@ subroutine stateInit(ini,phase,Nconstituents) volume - associate(stt => state(ins(phase))) + associate(stt => state(i)) if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume do e = 1,discretization_Nelems @@ -1794,21 +1777,22 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho(ph,me,ip,el) +pure function getRho(ph,me) - integer, intent(in) :: ph, me,ip,el + integer, intent(in) :: ph, me real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho + associate(prm => param(ins(ph))) - getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho = reshape(state(ins(ph))%rho(:,me),[prm%sum_N_sl,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) + ! 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%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & - getRho = 0.0_pReal + where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + getRho = 0.0_pReal end associate @@ -1819,21 +1803,22 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(ph,me,ip,el) +pure function getRho0(ph,me) - integer, intent(in) :: ph, me,ip,el + integer, intent(in) :: ph, me real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0 + associate(prm => param(ins(ph))) - getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(ins(ph))%rho(:,me),[prm%sum_N_sl,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) + ! 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%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & - getRho0 = 0.0_pReal + where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + getRho0 = 0.0_pReal end associate