fix: Use aperture table in poromechanics with conforming fractures (#3194)

paveltomin · Pavel Tomin · web-flow · commit 6c8202128980 · 2024-07-10T14:53:39.000-07:00
* decouple debug matrix output from logLevel

* Update SolidMechanicsLagrangeContact.cpp

* Update SinglePhasePoromechanicsConformingFractures.cpp

* Update SolidMechanicsLagrangeContact.cpp

* Update SurfaceElementSubRegion.hpp

* Update SinglePhaseBase.cpp

* Update HydrofractureSolver.cpp

* Update .integrated_tests.yaml

* Update BASELINE_NOTES.md

---------

Co-authored-by: Pavel Tomin &lt;“paveltomin@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-pr3006-5860-947a907
+  baseline: integratedTests/baseline_integratedTests-pr3194-6060-e5ba0ce
 
 allow_fail:
   all: ''
diff --git a/BASELINE_NOTES.md b/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.
diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
@@ -122,19 +122,20 @@ void ContactSolverBase::setFractureRegions( dataRepository::Group const & meshBo
 
 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 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
               break;
             }
           case FractureState::NewSlip:
+            {
+              newSlipCount += 1;
+              break;
+            }
           case FractureState::Slip:
             {
               slipCount += 1;
@@ -162,40 +167,43 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
     } );
 
     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 ) );
     } );
   }
 }
diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp b/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;
 
diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp
@@ -406,8 +406,8 @@ struct ComputeRotationMatricesKernel
     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];
diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
@@ -253,6 +253,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
 {
   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 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
         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 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
             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 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time,
 
   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 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
           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 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
                 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 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
   // 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,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/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,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/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();
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/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 >,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -750,12 +750,13 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
       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 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
           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 @@ void SinglePhaseBase::updateState( DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
 
-  if( m_keepFlowVariablesConstantDuringInitStep )
-    return;
-
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp
diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp
diff --git a/src/coreComponents/schema/docs/Hydrofracture.rst b/src/coreComponents/schema/docs/Hydrofracture.rst

Original file line number	Diff line number	Diff line change
`@@ -406,8 +406,8 @@ struct ComputeRotationMatricesKernel`
`406`	`406`	`forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )`
`407`	`407`	`{`
`408`	`408`
`409`		`- localIndex const & f0 = elemsToFaces[k][0];`
`410`		`- localIndex const & f1 = elemsToFaces[k][1];`
	`409`	`+ localIndex const f0 = elemsToFaces[k][0];`
	`410`	`+ localIndex const f1 = elemsToFaces[k][1];`
`411`	`411`
`412`	`412`	`real64 Nbar[3];`
`413`	`413`	`Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0];`
Original file line number	Diff line number	Diff line change
`@@ -253,6 +253,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain`
`253`	`253`	`{`
`254`	`254`	`GEOS_MARK_FUNCTION;`
`255`	`255`
	`256`	`+ real64 minNormalTractionTolerance( 1e10 );`
	`257`	`+ real64 maxNormalTractionTolerance( -1e10 );`
	`258`	`+ real64 minNormalDisplacementTolerance( 1e10 );`
	`259`	`+ real64 maxNormalDisplacementTolerance( -1e10 );`
	`260`	`+ real64 minSlidingTolerance( 1e10 );`
	`261`	`+ real64 maxSlidingTolerance( -1e10 );`
	`262`	`+`
`256`	`263`	`forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,`
`257`	`264`	`MeshLevel & mesh,`
`258`	`265`	`arrayView1d< string const > const & )`
`@@ -302,6 +309,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain`
`302`	`309`	`arrayView1d< real64 > const & slidingTolerance =`
`303`	`310`	`subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() );`
`304`	`311`
	`312`	`+ RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalTractionTolerance( 1e10 );`
	`313`	`+ RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalTractionTolerance( -1e10 );`
	`314`	`+ RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalDisplacementTolerance( 1e10 );`
	`315`	`+ RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalDisplacementTolerance( -1e10 );`
	`316`	`+ RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionSlidingTolerance( 1e10 );`
	`317`	`+ RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionSlidingTolerance( -1e10 );`
	`318`	`+`
`305`	`319`	`forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )`
`306`	`320`	`{`
`307`	`321`	`if( ghostRank[kfe] < 0 )`
`@@ -385,15 +399,36 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain`
`385`	`399`	`LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area );`
`386`	`400`
`387`	`401`	`// Finally, compute tolerances for the given fracture element`
	`402`	`+`
`388`	`403`	`normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7;`
	`404`	`+ minSubRegionNormalDisplacementTolerance.min( normalDisplacementTolerance[kfe] );`
	`405`	`+ maxSubRegionNormalDisplacementTolerance.max( normalDisplacementTolerance[kfe] );`
	`406`	`+`
`389`	`407`	`slidingTolerance[kfe] = sqrt( rotatedInvStiffApprox[ 1 ][ 1 ] * rotatedInvStiffApprox[ 1 ][ 1 ] +`
`390`	`408`	`rotatedInvStiffApprox[ 2 ][ 2 ] * rotatedInvStiffApprox[ 2 ][ 2 ] ) * averageYoungModulus / 2.e+7;`
	`409`	`+ minSubRegionSlidingTolerance.min( slidingTolerance[kfe] );`
	`410`	`+ maxSubRegionSlidingTolerance.max( slidingTolerance[kfe] );`
	`411`	`+`
`391`	`412`	`normalTractionTolerance[kfe] = 1.0 / 2.0 * averageConstrainedModulus / averageBoxSize0 * normalDisplacementTolerance[kfe];`
	`413`	`+ minSubRegionNormalTractionTolerance.min( normalTractionTolerance[kfe] );`
	`414`	`+ maxSubRegionNormalTractionTolerance.max( normalTractionTolerance[kfe] );`
`392`	`415`	`}`
`393`	`416`	`} );`
	`417`	`+`
	`418`	`+ minNormalDisplacementTolerance = std::min( minNormalDisplacementTolerance, minSubRegionNormalDisplacementTolerance.get() );`
	`419`	`+ maxNormalDisplacementTolerance = std::max( maxNormalDisplacementTolerance, maxSubRegionNormalDisplacementTolerance.get() );`
	`420`	`+ minSlidingTolerance = std::min( minSlidingTolerance, minSubRegionSlidingTolerance.get() );`
	`421`	`+ maxSlidingTolerance = std::max( maxSlidingTolerance, maxSubRegionSlidingTolerance.get() );`
	`422`	`+ minNormalTractionTolerance = std::min( minNormalTractionTolerance, minSubRegionNormalTractionTolerance.get() );`
	`423`	`+ maxNormalTractionTolerance = std::max( maxNormalTractionTolerance, maxSubRegionNormalTractionTolerance.get() );`
`394`	`424`	`}`
`395`	`425`	`} );`
`396`	`426`	`} );`
	`427`	`+`
	`428`	`+ GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{}: normal displacement tolerance = [{}, {}], sliding tolerance = [{}, {}], normal traction tolerance = [{}, {}]",`
	`429`	`+ this->getName(), minNormalDisplacementTolerance, maxNormalDisplacementTolerance,`
	`430`	`+ minSlidingTolerance, maxSlidingTolerance,`
	`431`	`+ minNormalTractionTolerance, maxNormalTractionTolerance ) );`
`397`	`432`	`}`
`398`	`433`
`399`	`434`	`void SolidMechanicsLagrangeContact::resetStateToBeginningOfStep( DomainPartition & domain )`
`@@ -586,7 +621,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time,`
`586`	`621`
`587`	`622`	`assembleContact( domain, dofManager, localMatrix, localRhs );`
`588`	`623`
`589`		`- // for sequenatial: add (fixed) pressure force contribution into residual (no derivatives)`
	`624`	`+ // for sequential: add (fixed) pressure force contribution into residual (no derivatives)`
`590`	`625`	`if( m_isFixedStressPoromechanicsUpdate )`
`591`	`626`	`{`
`592`	`627`	`forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,`
`@@ -2221,20 +2256,24 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai`
`2221`	`2256`	`if( ghostRank[kfe] < 0 )`
`2222`	`2257`	`{`
`2223`	`2258`	`integer const originalFractureState = fractureState[kfe];`
`2224`		`- if( originalFractureState == contact::FractureState::Open )`
	`2259`	`+ if( originalFractureState == FractureState::Open )`
`2225`	`2260`	`{`
`2226`		`- if( dispJump[kfe][0] > -normalDisplacementTolerance[kfe] )`
	`2261`	`+ if( dispJump[kfe][0] <= -normalDisplacementTolerance[kfe] )`
`2227`	`2262`	`{`
`2228`		`- fractureState[kfe] = contact::FractureState::Open;`
`2229`		`- }`
`2230`		`- else`
`2231`		`- {`
`2232`		`- fractureState[kfe] = contact::FractureState::Stick;`
	`2263`	`+ fractureState[kfe] = FractureState::Stick;`
	`2264`	`+ if( getLogLevel() >= 10 )`
	`2265`	`+ GEOS_LOG( GEOS_FMT( "{}: {} -> {}: dispJump = {}, normalDisplacementTolerance = {}",`
	`2266`	`+ kfe, originalFractureState, fractureState[kfe],`
	`2267`	`+ dispJump[kfe][0], normalDisplacementTolerance[kfe] ) );`
`2233`	`2268`	`}`
`2234`	`2269`	`}`
`2235`	`2270`	`else if( traction[kfe][0] > normalTractionTolerance[kfe] )`
`2236`	`2271`	`{`
`2237`		`- fractureState[kfe] = contact::FractureState::Open;`
	`2272`	`+ fractureState[kfe] = FractureState::Open;`
	`2273`	`+ if( getLogLevel() >= 10 )`
	`2274`	`+ GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}",`
	`2275`	`+ kfe, originalFractureState, fractureState[kfe],`
	`2276`	`+ traction[kfe][0], normalTractionTolerance[kfe] ) );`
`2238`	`2277`	`}`
`2239`	`2278`	`else`
`2240`	`2279`	`{`
`@@ -2245,29 +2284,33 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai`
`2245`	`2284`	`contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0],`
`2246`	`2285`	`dLimitTangentialTractionNorm_dTraction );`
`2247`	`2286`
`2248`		`- if( originalFractureState == contact::FractureState::Stick && currentTau >= limitTau )`
	`2287`	`+ if( originalFractureState == FractureState::Stick && currentTau >= limitTau )`
`2249`	`2288`	`{`
`2250`	`2289`	`currentTau *= (1.0 - m_slidingCheckTolerance);`
`2251`	`2290`	`}`
`2252`		`- else if( originalFractureState != contact::FractureState::Stick && currentTau <= limitTau )`
	`2291`	`+ else if( originalFractureState != FractureState::Stick && currentTau <= limitTau )`
`2253`	`2292`	`{`
`2254`	`2293`	`currentTau *= (1.0 + m_slidingCheckTolerance);`
`2255`	`2294`	`}`
`2256`	`2295`	`if( currentTau > limitTau )`
`2257`	`2296`	`{`
`2258`		`- if( originalFractureState == contact::FractureState::Stick )`
	`2297`	`+ if( originalFractureState == FractureState::Stick )`
`2259`	`2298`	`{`
`2260`		`- fractureState[kfe] = contact::FractureState::NewSlip;`
	`2299`	`+ fractureState[kfe] = FractureState::NewSlip;`
`2261`	`2300`	`}`
`2262`	`2301`	`else`
`2263`	`2302`	`{`
`2264`		`- fractureState[kfe] = contact::FractureState::Slip;`
	`2303`	`+ fractureState[kfe] = FractureState::Slip;`
`2265`	`2304`	`}`
`2266`	`2305`	`}`
`2267`	`2306`	`else`
`2268`	`2307`	`{`
`2269`		`- fractureState[kfe] = contact::FractureState::Stick;`
	`2308`	`+ fractureState[kfe] = FractureState::Stick;`
`2270`	`2309`	`}`
	`2310`	`+ if( getLogLevel() >= 10 && fractureState[kfe] != originalFractureState )`
	`2311`	`+ GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}",`
	`2312`	`+ kfe, originalFractureState, fractureState[kfe],`
	`2313`	`+ currentTau, limitTau ) );`
`2271`	`2314`	`}`
`2272`	`2315`
`2273`	`2316`	`changed += faceArea[kfe] * !compareFractureStates( originalFractureState, fractureState[kfe] );`
`@@ -2289,21 +2332,23 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai`
`2289`	`2332`	`// and total area of fracture elements`
`2290`	`2333`	`totalArea = MpiWrapper::sum( totalArea );`
`2291`	`2334`
	`2335`	`+ GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( " {}: changed area {} out of {}", getName(), changedArea, totalArea ) );`
	`2336`	`+`
`2292`	`2337`	`// Assume converged if changed area is below certain fraction of total area`
`2293`	`2338`	`return changedArea <= m_nonlinearSolverParameters.m_configurationTolerance * totalArea;`
`2294`	`2339`	`}`
`2295`	`2340`
`2296`	`2341`	`bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartition const & domain ) const`
`2297`	`2342`	`{`
`2298`		`- globalIndex numStick, numSlip, numOpen;`
	`2343`	`+ globalIndex numStick, numNewSlip, numSlip, numOpen;`
`2299`	`2344`	`forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,`
`2300`	`2345`	`MeshLevel const & mesh,`
`2301`	`2346`	`arrayView1d< string const > const & )`
`2302`	`2347`	`{`
`2303`		`- computeFractureStateStatistics( mesh, numStick, numSlip, numOpen );`
	`2348`	`+ computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen );`
`2304`	`2349`	`} );`
`2305`	`2350`
`2306`		`- return ( ( numSlip + numOpen ) == 0 );`
	`2351`	`+ return ( ( numNewSlip + numSlip + numOpen ) == 0 );`
`2307`	`2352`	`}`
`2308`	`2353`
`2309`	`2354`	`real64 SolidMechanicsLagrangeContact::setNextDt( real64 const & currentDt,`
Original file line number	Diff line number	Diff line change
`@@ -2585,9 +2585,6 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )`
`2585`	`2585`	`{`
`2586`	`2586`	`GEOS_MARK_FUNCTION;`
`2587`	`2587`
`2588`		`- if( m_keepFlowVariablesConstantDuringInitStep )`
`2589`		`- return;`
`2590`		`-`
`2591`	`2588`	`real64 maxDeltaPhaseVolFrac = 0.0;`
`2592`	`2589`	`forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,`
`2593`	`2590`	`MeshLevel & mesh,`
Original file line number	Diff line number	Diff line change
`@@ -750,12 +750,13 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,`
`750`	`750`	`singlePhaseBaseKernels::StatisticsKernel::`
`751`	`751`	`saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres );`
`752`	`752`
`753`		`- arrayView1d< real64 const > const dVol = subRegion.getField< fields::flow::deltaVolume >();`
	`753`	`+ arrayView1d< real64 > const dVol = subRegion.getField< fields::flow::deltaVolume >();`
`754`	`754`	`arrayView1d< real64 > const vol = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() );`
`755`	`755`
`756`	`756`	`forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )`
`757`	`757`	`{`
`758`	`758`	`vol[ei] += dVol[ei];`
	`759`	`+ dVol[ei] = 0.0;`
`759`	`760`	`} );`
`760`	`761`
`761`	`762`	`SingleFluidBase const & fluid =`
`@@ -803,7 +804,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,`
`803`	`804`	`if( volume[ei] * density_n[ei][0] > 1.1 * creationMass[ei] )`
`804`	`805`	`{`
`805`	`806`	`creationMass[ei] *= 0.75;`
`806`		`- if( creationMass[ei]<1.0e-20 )`
	`807`	`+ if( creationMass[ei] < 1.0e-20 )`
`807`	`808`	`{`
`808`	`809`	`creationMass[ei] = 0.0;`
`809`	`810`	`}`
`@@ -1255,9 +1256,6 @@ void SinglePhaseBase::updateState( DomainPartition & domain )`
`1255`	`1256`	`{`
`1256`	`1257`	`GEOS_MARK_FUNCTION;`
`1257`	`1258`
`1258`		`- if( m_keepFlowVariablesConstantDuringInitStep )`
`1259`		`- return;`
`1260`		`-`
`1261`	`1259`	`forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,`
`1262`	`1260`	`MeshLevel & mesh,`
`1263`	`1261`	`arrayView1d< string const > const & regionNames )`