feat: plastic strain output (#3384)

Author: ryar9534

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3486-9492-f0c817c
+  baseline: integratedTests/baseline_integratedTests-pr3384-9542-4e0f693
 allow_fail:
   all: ''
   streak: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3384 (2024-01-07)
+=====================
+Added plastic strain output.
+
 PR #3486 (2025-01-06)
 =====================
 useNewGravity became gravityDensityScheme.

src/coreComponents/constitutive/solid/ElasticIsotropic.hpp

@@ -135,6 +135,11 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates
                                  localIndex const q,
                                  real64 ( &elasticStrain )[6] ) const override final;
 
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrainInc( localIndex const k,
+                                    localIndex const q,
+                                    real64 ( &elasticStrainInc )[6] ) const override final;
+
   GEOS_HOST_DEVICE
   virtual real64 getBulkModulus( localIndex const k ) const override final
   {
@@ -171,6 +176,13 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates
 
 protected:
 
+  GEOS_HOST_DEVICE
+  virtual void computeElasticStrain( localIndex const k,
+                                     localIndex const q,
+                                     real64 const ( &stress )[6],
+                                     real64 ( &elasticStrainInc )[6] ) const;
+
+
   /// A reference to the ArrayView holding the bulk modulus for each element.
   arrayView1d< real64 const > const m_bulkModulus;
 
@@ -210,24 +222,53 @@ void ElasticIsotropicUpdates::getElasticStiffness( localIndex const k,
 }
 
 
+GEOS_HOST_DEVICE
+inline
+void ElasticIsotropicUpdates::computeElasticStrain( localIndex const k,
+                                                    localIndex const q,
+                                                    real64 const ( &stress )[6],
+                                                    real64 ( & elasticStrain)[6] ) const
+{
+  GEOS_UNUSED_VAR( q );
+  real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] );
+  real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] );
+
+  elasticStrain[0] = (    stress[0] - nu*stress[1] - nu*stress[2])/E;
+  elasticStrain[1] = (-nu*stress[0] +    stress[1] - nu*stress[2])/E;
+  elasticStrain[2] = (-nu*stress[0] - nu*stress[1] +    stress[2])/E;
+
+  elasticStrain[3] = stress[3] / m_shearModulus[k];
+  elasticStrain[4] = stress[4] / m_shearModulus[k];
+  elasticStrain[5] = stress[5] / m_shearModulus[k];
+}
+
 GEOS_HOST_DEVICE
 inline
 void ElasticIsotropicUpdates::getElasticStrain( localIndex const k,
                                                 localIndex const q,
                                                 real64 ( & elasticStrain)[6] ) const
 {
-  real64 const E = conversions::bulkModAndShearMod::toYoungMod( m_bulkModulus[k], m_shearModulus[k] );
-  real64 const nu = conversions::bulkModAndShearMod::toPoissonRatio( m_bulkModulus[k], m_shearModulus[k] );
 
-  elasticStrain[0] = (    m_newStress[k][q][0] - nu*m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E;
-  elasticStrain[1] = (-nu*m_newStress[k][q][0] +    m_newStress[k][q][1] - nu*m_newStress[k][q][2])/E;
-  elasticStrain[2] = (-nu*m_newStress[k][q][0] - nu*m_newStress[k][q][1] +    m_newStress[k][q][2])/E;
+  real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrain );
 
-  elasticStrain[3] = m_newStress[k][q][3] / m_shearModulus[k];
-  elasticStrain[4] = m_newStress[k][q][4] / m_shearModulus[k];
-  elasticStrain[5] = m_newStress[k][q][5] / m_shearModulus[k];
 }
 
+GEOS_HOST_DEVICE
+inline
+void ElasticIsotropicUpdates::getElasticStrainInc( localIndex const k,
+                                                   localIndex const q,
+                                                   real64 ( & elasticStrainInc)[6] ) const
+{
+
+  real64 stress[6] =
+  {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3],
+   m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrainInc );
+
+}
 
 GEOS_HOST_DEVICE
 inline

src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp

@@ -115,6 +115,11 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates
                                  localIndex const q,
                                  real64 ( &elasticStrain )[6] ) const override final;
 
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrainInc( localIndex const k,
+                                    localIndex const q,
+                                    real64 ( &elasticStrainInc )[6] ) const override final;
+
   GEOS_HOST_DEVICE
   virtual void viscousStateUpdate( localIndex const k,
                                    localIndex const q,
@@ -223,6 +228,60 @@ void ElasticIsotropicPressureDependentUpdates::getElasticStrain( localIndex cons
 }
 
 
+GEOS_HOST_DEVICE
+inline
+void ElasticIsotropicPressureDependentUpdates::getElasticStrainInc( localIndex const k,
+                                                                    localIndex const q,
+                                                                    real64 ( & elasticStrainInc)[6] ) const
+{
+  real64 const mu     = m_shearModulus[k];
+  real64 const p0     = m_refPressure;
+  real64 const eps_v0 = m_refStrainVol;
+  real64 const Cr     = m_recompressionIndex[k];
+  real64 deviator[6]{};
+  real64 stress[6]{};
+  real64 P;
+  real64 Q;
+
+  LvArray::tensorOps::copy< 6 >( stress, m_newStress[k][q] );
+
+  twoInvariant::stressDecomposition( stress,
+                                     P,
+                                     Q,
+                                     deviator );
+
+  real64 elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0;
+  real64 elasticStrainDev = Q/3./mu;
+
+  twoInvariant::strainRecomposition( elasticStrainVol,
+                                     elasticStrainDev,
+                                     deviator,
+                                     elasticStrainInc );
+
+  real64 oldStrain[6]{};
+
+  LvArray::tensorOps::copy< 6 >( stress, m_oldStress[k][q] );
+
+  twoInvariant::stressDecomposition( stress,
+                                     P,
+                                     Q,
+                                     deviator );
+
+  elasticStrainVol = LvArray::math::log( P/p0 ) * Cr * (-1.0) + eps_v0;
+  elasticStrainDev = Q/3./mu;
+
+  twoInvariant::strainRecomposition( elasticStrainVol,
+                                     elasticStrainDev,
+                                     deviator,
+                                     oldStrain );
+
+  for( localIndex i = 0; i<6; ++i )
+  {
+    elasticStrainInc[i] -= oldStrain[i];
+  }
+}
+
+
 GEOS_HOST_DEVICE
 inline
 void ElasticIsotropicPressureDependentUpdates::smallStrainUpdate( localIndex const k,

src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp

@@ -147,6 +147,16 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates
                                   real64 ( &stress )[6],
                                   DiscretizationOps & stiffness ) const final;
 
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrain( localIndex const k,
+                                 localIndex const q,
+                                 real64 ( &elasticStrain )[6] ) const override final;
+
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrainInc( localIndex const k,
+                                    localIndex const q,
+                                    real64 ( &elasticStrainInc )[6] ) const override final;
+
   // miscellaneous getters
 
   GEOS_HOST_DEVICE
@@ -162,6 +172,13 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates
     return LvArray::math::max( LvArray::math::max( m_c44[k], m_c55[k] ), m_c66[k] );
   }
 
+protected:
+  GEOS_HOST_DEVICE
+  virtual void computeElasticStrain( localIndex const k,
+                                     localIndex const q,
+                                     real64 const (&stress)[6],
+                                     real64 ( &elasticStrain )[6] ) const;
+
 private:
   /// A reference to the ArrayView holding c11 for each element.
   arrayView1d< real64 const > const m_c11;
@@ -219,6 +236,64 @@ void ElasticOrthotropicUpdates::getElasticStiffness( localIndex const k,
   stiffness[5][5] = m_c66[k];
 }
 
+
+GEOS_HOST_DEVICE
+inline
+void ElasticOrthotropicUpdates::computeElasticStrain( localIndex const k,
+                                                      localIndex const q,
+                                                      real64 const (&stress)[6],
+                                                      real64 ( & elasticStrain)[6] ) const
+{
+
+  GEOS_UNUSED_VAR( q );
+
+  real64 const detC = m_c11[k]*(m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k]) - m_c12[k]*(m_c12[k]*m_c33[k] - m_c23[k]*m_c13[k]) + m_c13[k]*(m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k]);
+
+  elasticStrain[0] =
+    ( (m_c22[k]*m_c33[k] - m_c23[k]*m_c23[k])*stress[0] + (m_c13[k]*m_c23[k] - m_c12[k]*m_c33[k])*stress[1] + (m_c12[k]*m_c23[k] - m_c13[k]*m_c22[k])*stress[2] ) /
+    detC;
+  elasticStrain[1] =
+    ( (m_c23[k]*m_c13[k] - m_c12[k]*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*m_c12[k] - m_c11[k]*m_c23[k])*stress[2] ) /
+    detC;
+  elasticStrain[2] =
+    ( (m_c12[k]*m_c23[k] - m_c22[k]*m_c13[k])*stress[0] + (m_c12[k]*m_c13[k] - m_c11[k]*m_c23[k])*stress[1] + (m_c11[k]*m_c22[k] - m_c12[k]*m_c12[k])*stress[2] ) /
+    detC;
+
+  elasticStrain[3] = stress[3] / m_c44[k];
+  elasticStrain[4] = stress[4] / m_c55[k];
+  elasticStrain[5] = stress[5] / m_c66[k];
+}
+
+
+GEOS_HOST_DEVICE
+inline
+void ElasticOrthotropicUpdates::getElasticStrain( localIndex const k,
+                                                  localIndex const q,
+                                                  real64 ( & elasticStrain)[6] ) const
+{
+
+  real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrain );
+
+}
+
+GEOS_HOST_DEVICE
+inline
+void ElasticOrthotropicUpdates::getElasticStrainInc( localIndex const k,
+                                                     localIndex const q,
+                                                     real64 ( & elasticStrainInc)[6] ) const
+{
+
+  real64 stress[6] =
+  {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3],
+   m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrainInc );
+
+}
+
+
 inline
 GEOS_HOST_DEVICE
 void ElasticOrthotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k,

src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp

@@ -144,6 +144,16 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates
   GEOS_HOST_DEVICE
   virtual void getElasticStiffness( localIndex const k, localIndex const q, real64 ( &stiffness )[6][6] ) const override final;
 
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrain( localIndex const k,
+                                 localIndex const q,
+                                 real64 ( &elasticStrain )[6] ) const override final;
+
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrainInc( localIndex const k,
+                                    localIndex const q,
+                                    real64 ( &elasticStrainInc )[6] ) const override final;
+
   /**
    * @brief Getter for apparent shear modulus.
    * @return reference to shear modulus that will be used for computing stabilization scalling parameter.
@@ -155,6 +165,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates
   }
 
 
+protected:
+
+  GEOS_HOST_DEVICE
+  virtual void computeElasticStrain( localIndex const k,
+                                     localIndex const q,
+                                     real64 const (&stress)[6],
+                                     real64 ( &elasticStrain )[6] ) const;
+
+
 private:
 
   /// A reference to the ArrayView holding c11 for each element.
@@ -200,6 +219,57 @@ void ElasticTransverseIsotropicUpdates::getElasticStiffness( localIndex const k,
   stiffness[5][5] = m_c66[k];
 }
 
+GEOS_HOST_DEVICE
+inline
+void ElasticTransverseIsotropicUpdates::computeElasticStrain( localIndex const k,
+                                                              localIndex const q,
+                                                              real64 const (&stress)[6],
+                                                              real64 ( & elasticStrain)[6] ) const
+{
+  GEOS_UNUSED_VAR( q );
+  real64 const c12 = ( m_c11[k] - 2.0 * m_c66[k] );
+  real64 const detC = m_c11[k]*(m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k]) - c12*(c12*m_c33[k] - m_c13[k]*m_c13[k]) + m_c13[k]*(c12*m_c13[k] - m_c11[k]*m_c13[k]);
+
+  elasticStrain[0] =
+    ( (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[0] + (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[1] + (c12*m_c13[k] - m_c13[k]*m_c11[k])*stress[2] ) / detC;
+  elasticStrain[1] =
+    ( (m_c13[k]*m_c13[k] - c12*m_c33[k])*stress[0] + (m_c11[k]*m_c33[k] - m_c13[k]*m_c13[k])*stress[1] + (m_c13[k]*c12 - m_c11[k]*m_c13[k])*stress[2] ) / detC;
+  elasticStrain[2] = ( (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[0] + (c12*m_c13[k] - m_c11[k]*m_c13[k])*stress[1] + (m_c11[k]*m_c11[k] - c12*c12)*stress[2] ) / detC;
+
+  elasticStrain[3] = stress[3] / m_c44[k];
+  elasticStrain[4] = stress[4] / m_c44[k];
+  elasticStrain[5] = stress[5] / m_c66[k];
+}
+
+
+GEOS_HOST_DEVICE
+inline
+void ElasticTransverseIsotropicUpdates::getElasticStrain( localIndex const k,
+                                                          localIndex const q,
+                                                          real64 ( & elasticStrain)[6] ) const
+{
+
+  real64 stress[6] = {m_newStress[k][q][0], m_newStress[k][q][1], m_newStress[k][q][2], m_newStress[k][q][3], m_newStress[k][q][4], m_newStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrain );
+
+}
+
+GEOS_HOST_DEVICE
+inline
+void ElasticTransverseIsotropicUpdates::getElasticStrainInc( localIndex const k,
+                                                             localIndex const q,
+                                                             real64 ( & elasticStrainInc)[6] ) const
+{
+
+  real64 stress[6] =
+  {m_newStress[k][q][0] - m_oldStress[k][q][0], m_newStress[k][q][1] - m_oldStress[k][q][1], m_newStress[k][q][2] - m_oldStress[k][q][2], m_newStress[k][q][3] - m_oldStress[k][q][3],
+   m_newStress[k][q][4] - m_oldStress[k][q][4], m_newStress[k][q][5] - m_oldStress[k][q][5]};
+
+  computeElasticStrain( k, q, stress, elasticStrainInc );
+
+}
+
 inline
 GEOS_HOST_DEVICE
 void ElasticTransverseIsotropicUpdates::smallStrainNoStateUpdate_StressOnly( localIndex const k,

src/coreComponents/constitutive/solid/SolidBase.hpp

@@ -337,6 +337,24 @@ class SolidBaseUpdates
     GEOS_ERROR( "getElasticStrain() not implemented for this model" );
   }
 
+  /**
+   * @brief Return the current elastic strain increment at a given material point (small-strain interface)
+   *
+   * @param k the element inex
+   * @param q the quadrature index
+   * @param elasticStrainInc Current elastic strain increment
+   */
+  GEOS_HOST_DEVICE
+  virtual void getElasticStrainInc( localIndex const k,
+                                    localIndex const q,
+                                    real64 ( & elasticStrainInc )[6] ) const
+  {
+    GEOS_UNUSED_VAR( k );
+    GEOS_UNUSED_VAR( q );
+    GEOS_UNUSED_VAR( elasticStrainInc );
+    GEOS_ERROR( "getElasticStrainInc() of SolidBase was called." );
+  }
+
   /**
    * @brief Perform a viscous (rate-dependent) state update
    *