feat: averaging of stress components using quadrature (#3279)

ryar9534 · Ryan Michael Aronson · Ryan Michael Aronson · web-flow · commit 6a90f8c5a8d3 · 2025-10-24T14:18:14.000-05:00
* initial attempt

* fixed compilation error, stresses match on simple problem

* renaming data layout for stress/strain and unifying a bit

* rename StrainHelper.hpp

* start collapsing stress and strain averaging into one kernel

* make the old stress data NOPLOT

* update SolidMechanicsFields

* add missing arguments

* change plasticStrain to averagePlasticStrain

* uncrustify

* make newStress NOPLOT

* update Kirsch problem to output and verify avg stress

* update python script accordingly

* update rst file for the definition of averageStress

* update .integrated_tests.yaml

* update BASELINE_NOTES.md

---------

Co-authored-by: Ryan Michael Aronson &lt;aronson1@ruby968.llnl.gov&gt;
Co-authored-by: Ryan Michael Aronson &lt;aronson1@ruby962.llnl.gov&gt;
Co-authored-by: Jian HUANG &lt;jian.huang@totalenergies.com&gt;
Co-authored-by: Jian Huang &lt;53012159+jhuang2601@users.noreply.github.com&gt;
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3851-14171-9d950d6
+  baseline: integratedTests/baseline_integratedTests-pr3279-14414-db55426
 
 allow_fail:
   all: ''
diff --git a/BASELINE_NOTES.md b/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 #3279 (2025-10-24) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3279-14414-db55426.tar.gz>
+=====================
+Output cell-wise average of each stress and strain component.
+
 PR #3851 (2025-10-13) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3851-14171-9d950d6.tar.gz>
 =====================
 Enable `FullyImplicit` for `SinglePhaseReservoirPoromechanicsConformingFractures`.
diff --git a/inputFiles/solidMechanics/KirschProblem_base.xml b/inputFiles/solidMechanics/KirschProblem_base.xml
@@ -96,7 +96,7 @@
     <PackCollection
       name="stressCollection"
       objectPath="ElementRegions/Omega/cb1"
-      fieldName="rock_stress"/>
+      fieldName="averageStress"/>
       
     <PackCollection
       name="displacementCollection"
@@ -112,11 +112,6 @@
       plotFileRoot="plot"
     />
 
-    <Silo
-      name="siloOutput"
-      plotLevel="3"
-    />
-
     <TimeHistory
       name="stressOutput"
       sources="{/Tasks/stressCollection}"
diff --git a/inputFiles/solidMechanics/KirschProblem_benchmark.xml b/inputFiles/solidMechanics/KirschProblem_benchmark.xml
@@ -40,11 +40,6 @@
       target="/Outputs/vtkOutput"
     />
 
-    <PeriodicEvent
-      name="outputs1"
-      target="/Outputs/siloOutput"
-    />
-
     <PeriodicEvent
       name="stressCollection"
       timeFrequency="1"
diff --git a/inputFiles/solidMechanics/KirschProblem_smoke.xml b/inputFiles/solidMechanics/KirschProblem_smoke.xml
@@ -39,11 +39,6 @@
       targetExactTimestep="1"
       target="/Outputs/vtkOutput"
     />
-
-    <PeriodicEvent
-      name="outputs1"
-      target="/Outputs/siloOutput"
-    />
     
     <PeriodicEvent
       name="restarts"
diff --git a/src/coreComponents/common/DataLayouts.hpp b/src/coreComponents/common/DataLayouts.hpp
@@ -155,24 +155,24 @@ namespace cells
 /// Cell node map permutation when using cuda.
 using NODE_MAP_PERMUTATION = RAJA::PERM_JI;
 
-/// Cell strain permutation when using cuda
-using STRAIN_PERM = RAJA::PERM_JI;
+/// Cell tensor (i.e. average stress and strain) permutation when using cuda
+using RANK2_TENSOR_PERM = RAJA::PERM_JI;
 
 #else
 
 /// Cell node map permutation when not using cuda.
 using NODE_MAP_PERMUTATION = RAJA::PERM_IJ;
 
-/// Cell strain permutation when not using cuda
-using STRAIN_PERM = RAJA::PERM_IJ;
+/// Cell tensor (i.e. average stress and strain) permutation when not using cuda
+using RANK2_TENSOR_PERM = RAJA::PERM_IJ;
 
 #endif
 
 /// Cell node map unit stride dimension.
 static constexpr int NODE_MAP_USD = LvArray::typeManipulation::getStrideOneDimension( NODE_MAP_PERMUTATION {} );
 
 /// Cell strain unit stride dimension
-static constexpr int STRAIN_USD = LvArray::typeManipulation::getStrideOneDimension( STRAIN_PERM {} );
+static constexpr int RANK2_TENSOR_USD = LvArray::typeManipulation::getStrideOneDimension( RANK2_TENSOR_PERM {} );
 
 } // namespace cells
 
diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp
@@ -38,7 +38,7 @@ DECLARE_FIELD( stress,
                "stress",
                array3dLayoutStress,
                0,
-               LEVEL_0,
+               NOPLOT,
                WRITE_AND_READ,
                "Current material stress" );
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsFields.hpp
@@ -42,9 +42,13 @@ using array2dLayoutIncrDisplacement = array2d< real64, nodes::INCR_DISPLACEMENT_
 using arrayView2dLayoutIncrDisplacement = arrayView2d< real64, nodes::INCR_DISPLACEMENT_USD >;
 using arrayViewConst2dLayoutIncrDisplacement = arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD >;
 
-using array2dLayoutStrain = array2d< real64, cells::STRAIN_PERM >;
-using arrayView2dLayoutStrain = arrayView2d< real64, cells::STRAIN_USD >;
-using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::STRAIN_USD >;
+using array2dLayoutStrain = array2d< real64, cells::RANK2_TENSOR_PERM >;
+using arrayView2dLayoutStrain = arrayView2d< real64, cells::RANK2_TENSOR_USD >;
+using arrayViewConst2dLayoutStrain = arrayView2d< real64 const, cells::RANK2_TENSOR_USD >;
+
+using array2dLayoutAvgStress = array2d< real64, cells::RANK2_TENSOR_PERM >;
+using arrayView2dLayoutAvgStress = arrayView2d< real64, cells::RANK2_TENSOR_USD >;
+using arrayViewConst2dLayoutAvgStress = arrayView2d< real64 const, cells::RANK2_TENSOR_USD >;
 
 using array2dLayoutVelocity = array2d< real64, nodes::VELOCITY_PERM >;
 using arrayView2dLayoutVelocity = arrayView2d< real64, nodes::VELOCITY_USD >;
@@ -71,16 +75,24 @@ DECLARE_FIELD( incrementalDisplacement,
                WRITE_AND_READ,
                "Incremental displacements for the current time step on the nodes" );
 
-DECLARE_FIELD( strain,
-               "strain",
+DECLARE_FIELD( averageStrain,
+               "averageStrain",
                array2dLayoutStrain,
                0,
                LEVEL_0,
                WRITE_AND_READ,
-               "Average strain in cell" );
+               "Quadrature averaged strain components in cell" );
+
+DECLARE_FIELD( averageStress,
+               "averageStress",
+               array2dLayoutAvgStress,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Quadrature averaged stress components in cell" );
 
-DECLARE_FIELD( plasticStrain,
-               "plasticStrain",
+DECLARE_FIELD( averagePlasticStrain,
+               "averagePlasticStrain",
                array2dLayoutStrain,
                0,
                LEVEL_0,
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
@@ -168,8 +168,9 @@ void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies )
                                                               [&]( localIndex const,
                                                                    ElementSubRegionBase & subRegion )
     {
-      subRegion.registerField< solidMechanics::strain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
-      subRegion.registerField< solidMechanics::plasticStrain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averageStrain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averagePlasticStrain >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
+      subRegion.registerField< solidMechanics::averageStress >( getName() ).setDimLabels( 1, voightLabels ).reference().resizeDimension< 1 >( 6 );
     } );
 
     NodeManager & nodes = meshLevel.getNodeManager();
@@ -924,29 +925,32 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS
       string const & solidMaterialName = subRegion.template getReference< string >( viewKeyStruct::solidMaterialNamesString() );
       SolidBase & constitutiveRelation = getConstitutiveModel< SolidBase >( subRegion, solidMaterialName );
 
+      arrayView3d< real64 const, solid::STRESS_USD > const stress = constitutiveRelation.getStress();
 
-      solidMechanics::arrayView2dLayoutStrain strain = subRegion.getField< solidMechanics::strain >();
-      solidMechanics::arrayView2dLayoutStrain plasticStrain = subRegion.getField< solidMechanics::plasticStrain >();
+      solidMechanics::arrayView2dLayoutStrain avgStrain = subRegion.getField< solidMechanics::averageStrain >();
+      solidMechanics::arrayView2dLayoutStrain avgPlasticStrain = subRegion.getField< solidMechanics::averagePlasticStrain >();
+      solidMechanics::arrayView2dLayoutAvgStress avgStress = subRegion.getField< solidMechanics::averageStress >();
 
       constitutive::ConstitutivePassThru< SolidBase >::execute( constitutiveRelation, [&] ( auto & solidModel )
       {
-
         using SOLID_TYPE = TYPEOFREF( solidModel );
 
         finiteElement::FiniteElementBase & subRegionFE = subRegion.template getReference< finiteElement::FiniteElementBase >( this->getDiscretizationName());
         finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::dispatch3D( subRegionFE, [&] ( auto const finiteElement )
         {
           using FE_TYPE = decltype( finiteElement );
-          AverageStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager,
-                                                                                                                          mesh.getEdgeManager(),
-                                                                                                                          mesh.getFaceManager(),
-                                                                                                                          subRegion,
-                                                                                                                          finiteElement,
-                                                                                                                          solidModel,
-                                                                                                                          disp,
-                                                                                                                          uhat,
-                                                                                                                          strain,
-                                                                                                                          plasticStrain );
+          AverageStressStrainOverQuadraturePointsKernelFactory::createAndLaunch< FE_TYPE, SOLID_TYPE, parallelDevicePolicy<> >( nodeManager,
+                                                                                                                                mesh.getEdgeManager(),
+                                                                                                                                mesh.getFaceManager(),
+                                                                                                                                subRegion,
+                                                                                                                                finiteElement,
+                                                                                                                                solidModel,
+                                                                                                                                disp,
+                                                                                                                                uhat,
+                                                                                                                                avgStrain,
+                                                                                                                                avgPlasticStrain,
+                                                                                                                                stress,
+                                                                                                                                avgStress );
         } );
 
 
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
@@ -23,7 +23,7 @@
 #include "common/format/EnumStrings.hpp"
 #include "common/TimingMacros.hpp"
 #include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
-#include "kernels/StrainHelper.hpp"
+#include "kernels/StressStrainAverageKernels.hpp"
 #include "mesh/mpiCommunications/CommunicationTools.hpp"
 #include "mesh/mpiCommunications/MPI_iCommData.hpp"
 #include "physicsSolvers/PhysicsSolverBase.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp
@@ -104,8 +104,8 @@ bool SolidMechanicsStateReset::execute( real64 const time_n,
                                                                 [&]( localIndex const,
                                                                      ElementSubRegionBase & subRegion )
       {
-        subRegion.getField< solidMechanics::strain >().zero();
-        subRegion.getField< solidMechanics::plasticStrain >().zero();
+        subRegion.getField< solidMechanics::averageStrain >().zero();
+        subRegion.getField< solidMechanics::averagePlasticStrain >().zero();
       } );
 
     }
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp
@@ -14,11 +14,11 @@
  */
 
 /**
- * @file StrainHelper.hpp
+ * @file StressStrainAverageKernels.hpp
  */
 
-#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_
-#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_
+#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_
+#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_
 
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
@@ -30,15 +30,17 @@
 
 namespace geos
 {
+
+
 /**
- * @class AverageStrainOverQuadraturePoints
+ * @class AverageStressStrainOverQuadraturePoints
  * @tparam SUBREGION_TYPE the subRegion type
  * @tparam FE_TYPE the finite element type
  * @tparam SOLID_TYPE the solid mechanics constitutuve type
  */
 template< typename FE_TYPE,
           typename SOLID_TYPE >
-class AverageStrainOverQuadraturePoints :
+class AverageStressStrainOverQuadraturePoints :
   public AverageOverQuadraturePointsBase< CellElementSubRegion,
                                           FE_TYPE >
 {
@@ -61,17 +63,21 @@ class AverageStrainOverQuadraturePoints :
    * @param finiteElementSpace the finite element space
    * @param displacement the displacement solution field
    * @param avgStrain the strain averaged over quadrature points
+   * @param stress the stress solution field
+   * @param avgStress the stress averaged over quadrature points
    */
-  AverageStrainOverQuadraturePoints( NodeManager & nodeManager,
-                                     EdgeManager const & edgeManager,
-                                     FaceManager const & faceManager,
-                                     CellElementSubRegion const & elementSubRegion,
-                                     FE_TYPE const & finiteElementSpace,
-                                     SOLID_TYPE const & solidModel,
-                                     fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
-                                     fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc,
-                                     fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
-                                     fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain ):
+  AverageStressStrainOverQuadraturePoints( NodeManager & nodeManager,
+                                           EdgeManager const & edgeManager,
+                                           FaceManager const & faceManager,
+                                           CellElementSubRegion const & elementSubRegion,
+                                           FE_TYPE const & finiteElementSpace,
+                                           SOLID_TYPE const & solidModel,
+                                           fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
+                                           fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc,
+                                           fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
+                                           fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain,
+                                           arrayView3d< real64 const, solid::STRESS_USD > const stress,
+                                           fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress ):
     Base( nodeManager,
           edgeManager,
           faceManager,
@@ -81,7 +87,9 @@ class AverageStrainOverQuadraturePoints :
     m_displacement( displacement ),
     m_displacementInc( displacementInc ),
     m_avgStrain( avgStrain ),
-    m_avgPlasticStrain( avgPlasticStrain )
+    m_avgPlasticStrain( avgPlasticStrain ),
+    m_stress( stress ),
+    m_avgStress( avgStress )
   {}
 
   /**
@@ -115,7 +123,7 @@ class AverageStrainOverQuadraturePoints :
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] = 0.0;
-      //m_avgPlasticStrain[k][icomp] = 0.0;
+      m_avgStress[k][icomp] = 0.0;
     }
   }
 
@@ -147,6 +155,7 @@ class AverageStrainOverQuadraturePoints :
     for( int icomp = 0; icomp < 6; ++icomp )
     {
       m_avgStrain[k][icomp] += conversionFactor[icomp]*detJxW*strain[icomp]/m_elementVolume[k];
+      m_avgStress[k][icomp] += detJxW*m_stress[k][q][icomp]/m_elementVolume[k];
 
       // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems
       // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless
@@ -201,15 +210,21 @@ class AverageStrainOverQuadraturePoints :
   /// The average plastic strain
   fields::solidMechanics::arrayView2dLayoutStrain const m_avgPlasticStrain;
 
+  /// The stress solution
+  arrayView3d< real64 const, solid::STRESS_USD > const m_stress;
+
+  /// The average stress
+  fields::solidMechanics::arrayView2dLayoutAvgStress const m_avgStress;
+
 };
 
 
 
 /**
- * @class AverageStrainOverQuadraturePointsKernelFactory
+ * @class AverageStressStrainOverQuadraturePointsKernelFactory
  * @brief Class to create and launch the kernel
  */
-class AverageStrainOverQuadraturePointsKernelFactory
+class AverageStressStrainOverQuadraturePointsKernelFactory
 {
 public:
 
@@ -240,13 +255,15 @@ class AverageStrainOverQuadraturePointsKernelFactory
                    fields::solidMechanics::arrayViewConst2dLayoutTotalDisplacement const displacement,
                    fields::solidMechanics::arrayViewConst2dLayoutIncrDisplacement const displacementInc,
                    fields::solidMechanics::arrayView2dLayoutStrain const avgStrain,
-                   fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain )
+                   fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain,
+                   arrayView3d< real64 const, solid::STRESS_USD > const stress,
+                   fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress )
   {
-    AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >
+    AverageStressStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >
     kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace,
-            solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain );
+            solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain, stress, avgStress );
 
-    AverageStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template
+    AverageStressStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template
     kernelLaunch< POLICY >( elementSubRegion.size(), kernel );
   }
 };
@@ -256,4 +273,5 @@ class AverageStrainOverQuadraturePointsKernelFactory
 }
 
 
-#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRAINHELPER_HPP_ */
+
+#endif /* GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_STRESSSTRAINAVERAGEKERNELS_HPP_ */
diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/Example.rst
@@ -140,6 +140,8 @@ In the ``Tasks`` section, ``PackCollection`` tasks are defined to collect time h
 Either the entire field or specified named sets of indices in the field can be collected. 
 In this example, ``stressCollection`` and ``displacementCollection`` tasks are specified to output the resultant stresses (tensor stored as an array with Voigt notation) and total displacement field (stored as a 3-component vector) respectively.
 
+Note: The output field named averageStress (comprising six components per cell) represents the cell-wise average of each stress component, calculated using the finite element quadrature rule. The same definition applies to averageStrain.
+
 .. literalinclude:: ../../../../../../../inputFiles/solidMechanics/KirschProblem_base.xml
     :language: xml
     :start-after: <!-- SPHINX_TASKS -->
diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/displacement_history.hdf5 b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/displacement_history.hdf5
diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/kirschWellboreFigure.py b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/kirschWellboreFigure.py
diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/stress_history.hdf5 b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/kirschWellbore/stress_history.hdf5
