feat: MGR Strategy for ALM - Face Bubble Functions for Wedge Elements…

… - Multipoint Integration Rules for Tetrahedra and Triangles (#3395) * Adding linear solving strategy * Adding bubble functions for wedge elements * Fixed face bubble functions for wedge elements * Adding parameters to tune the iterative penalty coefficients * Adding points for quadrature rule * Scaling dispJump by area * Added var to set the number of quadrature points for tetrahedron and triangle * Cleaning up SolidMechanicsALMSimultaneousKernels.hpp * Cleaning up H1_TriangleFace_Lagrange1_Gauss.hpp and H1_Tetrahedron_Lagrange1_Gauss.hpp * Fixing some review comments * Fixing compilation issues * Fixed bug (int-real type casting) * Modified the default input parameters
GEOS-DEV · Jan 24, 2025 · 79c0422 · 79c0422
1 parent 50504b5
commit 79c0422
Show file tree

Hide file tree

Showing 32 changed files with 863 additions and 189 deletions.
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3416-9790-5bdb1fa
+  baseline: integratedTests/baseline_integratedTests-pr3395-9832-31f145e
 
 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  #3395 (2024-01-22)
+=====================
+Add new fields and change the default input for some tests.
+
 PR #3416 (2024-01-21)
 =====================
 Refactoring of induced seismicity EQ solvers to add coupling.

diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml
@@ -64,7 +64,7 @@
     <SurfaceGenerator
       name="SurfaceGen"
       targetRegions="{ Domain }"
-      rockToughness="1.0"
+      initialRockToughness="1.0"
       mpiCommOrder="1"
       fractureRegion="Fault"/>  
 
@@ -136,4 +136,4 @@
       targetExactTimestep="0"
       target="/Outputs/timeHistoryOutput"/>
   </Events>
-</Problem>
+</Problem>
diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml
@@ -65,7 +65,7 @@
     <SurfaceGenerator
       name="SurfaceGen"
       targetRegions="{ Domain }"
-      rockToughness="1.0"
+      initialRockToughness="1.0"
       mpiCommOrder="1"
       fractureRegion="Fault"/>  
 
@@ -137,4 +137,4 @@
       targetExactTimestep="0"
       target="/Outputs/timeHistoryOutput"/>
   </Events>
-</Problem>
+</Problem>
diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml
@@ -45,7 +45,7 @@
     <SurfaceGenerator
       name="SurfaceGen"
       targetRegions="{ Domain }"
-      rockToughness="1.0"
+      initialRockToughness="1.0"
       mpiCommOrder="1"
       fractureRegion="Fault"/>  
   </Solvers>
@@ -248,4 +248,4 @@
       targetExactTimestep="0"
       target="/Outputs/restart"/>   -->
   </Events> 
-</Problem>
+</Problem>
diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp
@@ -103,7 +103,6 @@ class CoulombFrictionUpdates : public FrictionBaseUpdates
                                arraySlice1d< real64 const > const & dispJump,
                                arraySlice1d< real64 const > const & penalty,
                                arraySlice1d< real64 const > const & traction,
-                               real64 const faceArea,
                                bool const symmetric,
                                bool const fixedLimitTau,
                                real64 const normalTractionTolerance,
@@ -119,7 +118,6 @@ class CoulombFrictionUpdates : public FrictionBaseUpdates
                                    arraySlice1d< real64 const > const & deltaDispJump,
                                    arraySlice1d< real64 const > const & penalty,
                                    arraySlice1d< real64 const > const & traction,
-                                   real64 const faceArea,
                                    arraySlice1d< real64 > const & tractionNew ) const override final;
 
   GEOS_HOST_DEVICE
@@ -363,7 +361,6 @@ inline void CoulombFrictionUpdates::updateTraction( arraySlice1d< real64 const >
                                                     arraySlice1d< real64 const > const & dispJump,
                                                     arraySlice1d< real64 const > const & penalty,
                                                     arraySlice1d< real64 const > const & traction,
-                                                    real64 const faceArea,
                                                     bool const symmetric,
                                                     bool const fixedLimitTau,
                                                     real64 const normalTractionTolerance,
@@ -380,9 +377,9 @@ inline void CoulombFrictionUpdates::updateTraction( arraySlice1d< real64 const >
 
   // Compute the trial traction
   real64 tractionTrial[ 3 ];
-  tractionTrial[ 0 ] = traction[0] + penalty[0] * dispJump[0] * faceArea;
-  tractionTrial[ 1 ] = traction[1] + penalty[1] * (dispJump[1] - oldDispJump[1]) * faceArea;
-  tractionTrial[ 2 ] = traction[2] + penalty[1] * (dispJump[2] - oldDispJump[2]) * faceArea;
+  tractionTrial[ 0 ] = traction[0] + penalty[0] * dispJump[0];
+  tractionTrial[ 1 ] = traction[1] + penalty[1] * (dispJump[1] - oldDispJump[1]);
+  tractionTrial[ 2 ] = traction[2] + penalty[1] * (dispJump[2] - oldDispJump[2]);
 
   // Compute tangential trial traction norm
   real64 const tau[2] = { tractionTrial[1],
@@ -502,16 +499,15 @@ inline void CoulombFrictionUpdates::updateTractionOnly( arraySlice1d< real64 con
                                                         arraySlice1d< real64 const > const & deltaDispJump,
                                                         arraySlice1d< real64 const > const & penalty,
                                                         arraySlice1d< real64 const > const & traction,
-                                                        real64 const faceArea,
                                                         arraySlice1d< real64 > const & tractionNew ) const
 {
 
   // TODO: Pass this tol as an argument or define a new class member
   real64 const zero = LvArray::NumericLimits< real64 >::epsilon;
 
-  tractionNew[0] = traction[0] + penalty[0] * dispJump[0] * faceArea;
-  tractionNew[1] = traction[1] + penalty[1] * deltaDispJump[1] * faceArea;
-  tractionNew[2] = traction[2] + penalty[1] * deltaDispJump[2] * faceArea;
+  tractionNew[0] = traction[0] + penalty[0] * dispJump[0];
+  tractionNew[1] = traction[1] + penalty[1] * deltaDispJump[1];
+  tractionNew[2] = traction[2] + penalty[1] * deltaDispJump[2];
 
   real64 const tau[2] = { tractionNew[1],
                           tractionNew[2] };

diff --git a/src/coreComponents/constitutive/contact/FrictionBase.hpp b/src/coreComponents/constitutive/contact/FrictionBase.hpp
@@ -125,7 +125,6 @@ class FrictionBaseUpdates
                                arraySlice1d< real64 const > const & dispJump,
                                arraySlice1d< real64 const > const & penalty,
                                arraySlice1d< real64 const > const & traction,
-                               real64 const faceArea,
                                bool const symmetric,
                                bool const fixedLimitTau,
                                real64 const normalTractionTolerance,
@@ -134,7 +133,7 @@ class FrictionBaseUpdates
                                real64 ( & tractionNew )[3],
                                integer & fractureState ) const
   {
-    GEOS_UNUSED_VAR( oldDispJump, dispJump, penalty, traction, faceArea, symmetric, fixedLimitTau,
+    GEOS_UNUSED_VAR( oldDispJump, dispJump, penalty, traction, symmetric, fixedLimitTau,
                      normalTractionTolerance, tangentialTractionTolerance,
                      dTraction_dDispJump, tractionNew, fractureState );
   }
@@ -153,9 +152,8 @@ class FrictionBaseUpdates
                                    arraySlice1d< real64 const > const & deltaDispJump,
                                    arraySlice1d< real64 const > const & penalty,
                                    arraySlice1d< real64 const > const & traction,
-                                   real64 const faceArea,
                                    arraySlice1d< real64 > const & tractionNew ) const
-  { GEOS_UNUSED_VAR( dispJump, deltaDispJump, penalty, traction, faceArea, tractionNew ); }
+  { GEOS_UNUSED_VAR( dispJump, deltaDispJump, penalty, traction, tractionNew ); }
 
   /**
    * @brief Check for the constraint satisfaction

diff --git a/src/coreComponents/finiteElement/CMakeLists.txt b/src/coreComponents/finiteElement/CMakeLists.txt
@@ -32,8 +32,8 @@ set( finiteElement_headers
      elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp
      elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp
      elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp
-     elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp
-     elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp
+     elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp
+     elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp
      elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp
      elementFormulations/ConformingVirtualElementOrder1.hpp
      elementFormulations/ConformingVirtualElementOrder1_impl.hpp

diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp b/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp
@@ -49,6 +49,11 @@ FiniteElementDiscretization::FiniteElementDiscretization( string const & name, G
     setInputFlag( InputFlags::OPTIONAL ).
     setApplyDefaultValue( 0 ).
     setDescription( "Specifier to indicate whether to force the use of VEM" );
+
+  registerWrapper( viewKeyStruct::useHighOrderQuadratureRuleString(), &m_useHighOrderQuadratureRule ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( 0 ).
+    setDescription( "Specifier to indicate whether to use a high order quadrature rule" );
 }
 
 FiniteElementDiscretization::~FiniteElementDiscretization()
@@ -68,7 +73,15 @@ FiniteElementDiscretization::factory( ElementType const parentElementShape ) con
   {
     switch( parentElementShape )
     {
-      case ElementType::Triangle:      return std::make_unique< H1_TriangleFace_Lagrange1_Gauss1 >();
+      case ElementType::Triangle:
+        if( m_useHighOrderQuadratureRule == 1 )
+        {
+          return std::make_unique< H1_TriangleFace_Lagrange1_Gauss4 >();
+        }
+        else
+        {
+          return std::make_unique< H1_TriangleFace_Lagrange1_Gauss1 >();
+        }
       case ElementType::Quadrilateral: return std::make_unique< H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >();
       // On polyhedra where FEM are available, we use VEM only if useVirtualElements is set to 1 in
       // the input file.
@@ -80,7 +93,14 @@ FiniteElementDiscretization::factory( ElementType const parentElementShape ) con
         }
         else
         {
-          return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss1 >();
+          if( m_useHighOrderQuadratureRule == 1 )
+          {
+            return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss14 >();
+          }
+          else
+          {
+            return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss1 >();
+          }
         }
       }
       case ElementType::Pyramid:

diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp b/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp
@@ -97,6 +97,7 @@ class FiniteElementDiscretization : public dataRepository::Group
     static constexpr char const * orderString() { return "order"; }
     static constexpr char const * formulationString() { return "formulation"; }
     static constexpr char const * useVemString() { return "useVirtualElements"; }
+    static constexpr char const * useHighOrderQuadratureRuleString() { return "useHighOrderQuadratureRule"; }
   };
 
   /// The order of the finite element basis
@@ -108,6 +109,9 @@ class FiniteElementDiscretization : public dataRepository::Group
   /// Optional parameter indicating if the class should use Virtual Elements.
   int m_useVem;
 
+  /// Optional parameter indicating if the class should use a high order quadrature rule.
+  int m_useHighOrderQuadratureRule;
+
   void postInputInitialization() override final;
 
 };

diff --git a/src/coreComponents/finiteElement/FiniteElementDispatch.hpp b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp
@@ -24,19 +24,20 @@
 #include "elementFormulations/ConformingVirtualElementOrder1.hpp"
 #include "elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp"
 #include "elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp"
-#include "elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp"
+#include "elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp"
 #include "elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp"
 #if !defined( GEOS_USE_HIP )
 #include "elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp"
 #endif
 #include "elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp"
-#include "elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp"
+#include "elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp"
 #include "LvArray/src/system.hpp"
 
 #define FE_1_TYPES \
   finiteElement::H1_Hexahedron_Lagrange1_GaussLegendre2, \
   finiteElement::H1_Wedge_Lagrange1_Gauss6, \
   finiteElement::H1_Tetrahedron_Lagrange1_Gauss1, \
+  finiteElement::H1_Tetrahedron_Lagrange1_Gauss14, \
   finiteElement::H1_Pyramid_Lagrange1_Gauss5
 
 #define GL_FE_TYPES \
@@ -88,7 +89,8 @@
 
 #define FE_TYPES_2D \
   finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2, \
-  finiteElement::H1_TriangleFace_Lagrange1_Gauss1
+  finiteElement::H1_TriangleFace_Lagrange1_Gauss1, \
+  finiteElement::H1_TriangleFace_Lagrange1_Gauss4
 
 #define BASE_FE_TYPES_2D FE_TYPES_2D