preparing data handling for explicit forward of flux

src/constitutive_plastic_nonlocal.f90

@@ -180,7 +180,8 @@ submodule(constitutive) plastic_nonlocal
   type(tNonlocalState), allocatable, dimension(:) :: &
     deltaState, &
     dotState, &
-    state
+    state, &
+    state0
 
   type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstance)
 
@@ -238,6 +239,7 @@ 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))
@@ -482,6 +484,7 @@ module subroutine plastic_nonlocal_init
 
     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,:)
@@ -1713,7 +1716,6 @@ 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 &
@@ -1931,6 +1933,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
 !--------------------------------------------------------------------------------------------------