fix: Use aperture table in poromechanics with conforming fractures

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3006-5860-947a907
+  baseline: integratedTests/baseline_integratedTests-pr3194-6060-e5ba0ce
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,18 +6,26 @@ 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 #3194 (2024-07-10)
+======================
+Use aperture table in poromechanics with conforming fractures. Rebaseline the corresponding cases.
+
+
 PR #3006 (2024-07-01)
 ======================
 Added baselines for new tests. Relaxing tolerances for singlePhasePoromechanics_FaultModel_smoke.
 
+
 PR #3196 (2024-06-28)
 ======================
 Added isLaggingFractureStencilWeightsUpdate to hydrofracture solve. Rebaseline because of the new input.
 
+
 PR #3177 (2024-06-28)
 ======================
 Added logLevel to TimeHistoryOutput. Rebaseline because of the new input flag.
 
+
 PR #3181 (2024-06-25)
 ======================
 Decouple debug matrix output from logLevel. Rebaseline because of the new input flag.

src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp

@@ -122,19 +122,20 @@
 
 void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
                                                         globalIndex & numStick,
+                                                        globalIndex & numNewSlip,
                                                         globalIndex & numSlip,
                                                         globalIndex & numOpen ) const
 {
   ElementRegionManager const & elemManager = mesh.getElemManager();
 
-  array1d< globalIndex > localCounter( 3 );
+  array1d< globalIndex > localCounter( 4 );
 
   elemManager.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion const & subRegion )
   {
     arrayView1d< integer const > const & ghostRank = subRegion.ghostRank();
     arrayView1d< integer const > const & fractureState = subRegion.getField< fields::contact::fractureState >();
 
-    RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), slipCount( 0 ), openCount( 0 );
+    RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), newSlipCount( 0 ), slipCount( 0 ), openCount( 0 );
     forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )
     {
       if( ghostRank[kfe] < 0 )
@@ -147,6 +148,10 @@
               break;
             }
           case FractureState::NewSlip:
+            {
+              newSlipCount += 1;
+              break;
+            }
           case FractureState::Slip:
             {
               slipCount += 1;
@@ -162,40 +167,43 @@
     } );
 
     localCounter[0] += stickCount.get();
-    localCounter[1] += slipCount.get();
-    localCounter[2] += openCount.get();
+    localCounter[1] += newSlipCount.get();
+    localCounter[2] += slipCount.get();
+    localCounter[3] += openCount.get();
   } );
 
-  array1d< globalIndex > totalCounter( 3 );
+  array1d< globalIndex > totalCounter( 4 );
 
   MpiWrapper::allReduce( localCounter.data(),
                          totalCounter.data(),
-                         3,
+                         4,
                          MPI_SUM,
                          MPI_COMM_GEOSX );
 
-  numStick = totalCounter[0];
-  numSlip  = totalCounter[1];
-  numOpen  = totalCounter[2];
+  numStick    = totalCounter[0];
+  numNewSlip  = totalCounter[1];
+  numSlip     = totalCounter[2];
+  numOpen     = totalCounter[3];
 }
 
 void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & domain ) const
 {
   if( getLogLevel() >=1 )
   {
     globalIndex numStick = 0;
+    globalIndex numNewSlip  = 0;
     globalIndex numSlip  = 0;
     globalIndex numOpen  = 0;
 
     forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                  MeshLevel const & mesh,
                                                                  arrayView1d< string const > const & )
     {
-      computeFractureStateStatistics( mesh, numStick, numSlip, numOpen );
+      computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen );
 
       GEOS_LOG_RANK_0( GEOS_FMT( "  Number of element for each fracture state:"
-                                 " stick: {:12} | slip:  {:12} | open:  {:12}",
-                                 numStick, numSlip, numOpen ) );
+                                 " stick: {:12} | new slip: {:12} | slip:  {:12} | open:  {:12}",
+                                 numStick, numNewSlip, numSlip, numOpen ) );
     } );
   }
 }

src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp

@@ -62,6 +62,7 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM
 
   void computeFractureStateStatistics( MeshLevel const & mesh,
                                        globalIndex & numStick,
+                                       globalIndex & numNewSlip,
                                        globalIndex & numSlip,
                                        globalIndex & numOpen ) const;
 

src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp

@@ -406,8 +406,8 @@
     forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
     {
 
-      localIndex const & f0 = elemsToFaces[k][0];
-      localIndex const & f1 = elemsToFaces[k][1];
+      localIndex const f0 = elemsToFaces[k][0];
+      localIndex const f1 = elemsToFaces[k][1];
 
       real64 Nbar[3];
       Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0];

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp

@@ -253,6 +253,13 @@
 {
   GEOS_MARK_FUNCTION;
 
+  real64 minNormalTractionTolerance( 1e10 );
+  real64 maxNormalTractionTolerance( -1e10 );
+  real64 minNormalDisplacementTolerance( 1e10 );
+  real64 maxNormalDisplacementTolerance( -1e10 );
+  real64 minSlidingTolerance( 1e10 );
+  real64 maxSlidingTolerance( -1e10 );
+
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel & mesh,
                                                                 arrayView1d< string const > const & )
@@ -302,6 +309,13 @@
         arrayView1d< real64 > const & slidingTolerance =
           subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() );
 
+        RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalTractionTolerance( 1e10 );
+        RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalTractionTolerance( -1e10 );
+        RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalDisplacementTolerance( 1e10 );
+        RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalDisplacementTolerance( -1e10 );
+        RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionSlidingTolerance( 1e10 );
+        RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionSlidingTolerance( -1e10 );
+
         forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )
         {
           if( ghostRank[kfe] < 0 )
@@ -385,15 +399,36 @@
             LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area );
 
             // Finally, compute tolerances for the given fracture element
+
             normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7;
+            minSubRegionNormalDisplacementTolerance.min( normalDisplacementTolerance[kfe] );
+            maxSubRegionNormalDisplacementTolerance.max( normalDisplacementTolerance[kfe] );
+
             slidingTolerance[kfe] = sqrt( rotatedInvStiffApprox[ 1 ][ 1 ] * rotatedInvStiffApprox[ 1 ][ 1 ] +
                                           rotatedInvStiffApprox[ 2 ][ 2 ] * rotatedInvStiffApprox[ 2 ][ 2 ] ) * averageYoungModulus / 2.e+7;
+            minSubRegionSlidingTolerance.min( slidingTolerance[kfe] );
+            maxSubRegionSlidingTolerance.max( slidingTolerance[kfe] );
+
             normalTractionTolerance[kfe] = 1.0 / 2.0 * averageConstrainedModulus / averageBoxSize0 * normalDisplacementTolerance[kfe];
+            minSubRegionNormalTractionTolerance.min( normalTractionTolerance[kfe] );
+            maxSubRegionNormalTractionTolerance.max( normalTractionTolerance[kfe] );
           }
         } );
+
+        minNormalDisplacementTolerance = std::min( minNormalDisplacementTolerance, minSubRegionNormalDisplacementTolerance.get() );
+        maxNormalDisplacementTolerance = std::max( maxNormalDisplacementTolerance, maxSubRegionNormalDisplacementTolerance.get() );
+        minSlidingTolerance = std::min( minSlidingTolerance, minSubRegionSlidingTolerance.get() );
+        maxSlidingTolerance = std::max( maxSlidingTolerance, maxSubRegionSlidingTolerance.get() );
+        minNormalTractionTolerance = std::min( minNormalTractionTolerance, minSubRegionNormalTractionTolerance.get() );
+        maxNormalTractionTolerance = std::max( maxNormalTractionTolerance, maxSubRegionNormalTractionTolerance.get() );
       }
     } );
   } );
+
+  GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{}: normal displacement tolerance = [{}, {}], sliding tolerance = [{}, {}], normal traction tolerance = [{}, {}]",
+                                      this->getName(), minNormalDisplacementTolerance, maxNormalDisplacementTolerance,
+                                      minSlidingTolerance, maxSlidingTolerance,
+                                      minNormalTractionTolerance, maxNormalTractionTolerance ) );
 }
 
 void SolidMechanicsLagrangeContact::resetStateToBeginningOfStep( DomainPartition & domain )
@@ -586,7 +621,7 @@
 
   assembleContact( domain, dofManager, localMatrix, localRhs );
 
-  // for sequenatial: add (fixed) pressure force contribution into residual (no derivatives)
+  // for sequential: add (fixed) pressure force contribution into residual (no derivatives)
   if( m_isFixedStressPoromechanicsUpdate )
   {
     forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
@@ -2221,20 +2256,24 @@
           if( ghostRank[kfe] < 0 )
           {
             integer const originalFractureState = fractureState[kfe];
-            if( originalFractureState == contact::FractureState::Open )
+            if( originalFractureState == FractureState::Open )
             {
-              if( dispJump[kfe][0] > -normalDisplacementTolerance[kfe] )
+              if( dispJump[kfe][0] <= -normalDisplacementTolerance[kfe] )
               {
-                fractureState[kfe] = contact::FractureState::Open;
-              }
-              else
-              {
-                fractureState[kfe] = contact::FractureState::Stick;
+                fractureState[kfe] = FractureState::Stick;
+                if( getLogLevel() >= 10 )
+                  GEOS_LOG( GEOS_FMT( "{}: {} -> {}: dispJump = {}, normalDisplacementTolerance = {}",
+                                      kfe, originalFractureState, fractureState[kfe],
+                                      dispJump[kfe][0], normalDisplacementTolerance[kfe] ) );
               }
             }
             else if( traction[kfe][0] > normalTractionTolerance[kfe] )
             {
-              fractureState[kfe] = contact::FractureState::Open;
+              fractureState[kfe] = FractureState::Open;
+              if( getLogLevel() >= 10 )
+                GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}",
+                                    kfe, originalFractureState, fractureState[kfe],
+                                    traction[kfe][0], normalTractionTolerance[kfe] ) );
             }
             else
             {
@@ -2245,29 +2284,33 @@
                 contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0],
                                                                    dLimitTangentialTractionNorm_dTraction );
 
-              if( originalFractureState == contact::FractureState::Stick && currentTau >= limitTau )
+              if( originalFractureState == FractureState::Stick && currentTau >= limitTau )
               {
                 currentTau *= (1.0 - m_slidingCheckTolerance);
               }
-              else if( originalFractureState != contact::FractureState::Stick && currentTau <= limitTau )
+              else if( originalFractureState != FractureState::Stick && currentTau <= limitTau )
               {
                 currentTau *= (1.0 + m_slidingCheckTolerance);
               }
               if( currentTau > limitTau )
               {
-                if( originalFractureState == contact::FractureState::Stick )
+                if( originalFractureState == FractureState::Stick )
                 {
-                  fractureState[kfe] = contact::FractureState::NewSlip;
+                  fractureState[kfe] = FractureState::NewSlip;
                 }
                 else
                 {
-                  fractureState[kfe] = contact::FractureState::Slip;
+                  fractureState[kfe] = FractureState::Slip;
                 }
               }
               else
               {
-                fractureState[kfe] = contact::FractureState::Stick;
+                fractureState[kfe] = FractureState::Stick;
               }
+              if( getLogLevel() >= 10 && fractureState[kfe] != originalFractureState )
+                GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}",
+                                    kfe, originalFractureState, fractureState[kfe],
+                                    currentTau, limitTau ) );
             }
 
             changed += faceArea[kfe] * !compareFractureStates( originalFractureState, fractureState[kfe] );
@@ -2289,21 +2332,23 @@
   // and total area of fracture elements
   totalArea = MpiWrapper::sum( totalArea );
 
+  GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "  {}: changed area {} out of {}", getName(), changedArea, totalArea ) );
+
   // Assume converged if changed area is below certain fraction of total area
   return changedArea <= m_nonlinearSolverParameters.m_configurationTolerance * totalArea;
 }
 
 bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartition const & domain ) const
 {
-  globalIndex numStick, numSlip, numOpen;
+  globalIndex numStick, numNewSlip, numSlip, numOpen;
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel const & mesh,
                                                                 arrayView1d< string const > const & )
   {
-    computeFractureStateStatistics( mesh, numStick, numSlip, numOpen );
+    computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen );
   } );
 
-  return ( ( numSlip + numOpen ) == 0 );
+  return ( ( numNewSlip + numSlip + numOpen ) == 0 );
 }
 
 real64 SolidMechanicsLagrangeContact::setNextDt( real64 const & currentDt,

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

@@ -2585,9 +2585,6 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
 
-  if( m_keepFlowVariablesConstantDuringInitStep )
-    return;
-
   real64 maxDeltaPhaseVolFrac = 0.0;
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

@@ -185,9 +185,6 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies )
       subRegion.registerField< fields::flow::hydraulicAperture >( getName() ).
         setApplyDefaultValue( faceRegion.getDefaultAperture() );
 
-      subRegion.registerField< fields::flow::minimumHydraulicAperture >( getName() ).
-        setApplyDefaultValue( faceRegion.getDefaultAperture() );
-
     } );
 
     FaceManager & faceManager = mesh.getFaceManager();

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp

@@ -168,14 +168,6 @@ DECLARE_FIELD( hydraulicAperture,
                WRITE_AND_READ,
                "Hydraulic aperture" );
 
-DECLARE_FIELD( minimumHydraulicAperture,
-               "minimumHydraulicAperture",
-               array1d< real64 >,
-               0,
-               LEVEL_0,
-               WRITE_AND_READ,
-               "minimum value of the hydraulic aperture" );
-
 DECLARE_FIELD( gravityCoefficient,
                "gravityCoefficient",
                array1d< real64 >,

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp

@@ -750,12 +750,13 @@
       singlePhaseBaseKernels::StatisticsKernel::
         saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres );
 
-      arrayView1d< real64 const > const dVol = subRegion.getField< fields::flow::deltaVolume >();
+      arrayView1d< real64 > const dVol = subRegion.getField< fields::flow::deltaVolume >();
       arrayView1d< real64 > const vol = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() );
 
       forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
       {
         vol[ei] += dVol[ei];
+        dVol[ei] = 0.0;
       } );
 
       SingleFluidBase const & fluid =
@@ -803,7 +804,7 @@
           if( volume[ei] * density_n[ei][0] > 1.1 * creationMass[ei] )
           {
             creationMass[ei] *= 0.75;
-            if( creationMass[ei]<1.0e-20 )
+            if( creationMass[ei] < 1.0e-20 )
             {
               creationMass[ei] = 0.0;
             }
@@ -1255,9 +1256,6 @@
 {
   GEOS_MARK_FUNCTION;
 
-  if( m_keepFlowVariablesConstantDuringInitStep )
-    return;
-
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )