fix: delta volume fix plus some other stuff

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3339-8707-7c55c70
+  baseline: integratedTests/baseline_integratedTests-pr3439-8742-0bffafb
 allow_fail:
   all: ''
   streak: ''

BASELINE_NOTES.md

@@ -6,6 +6,9 @@ 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 #3439 (2024-11-14)
+=====================
+EDFM bugfixes: derivatives sign, frac/cell element volume.
 
 PR #3339 (2024-11-14)
 =====================

inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml

@@ -13,7 +13,9 @@
       <NonlinearSolverParameters
         newtonTol="1.0e-1"
         newtonMaxIter="10"
-        maxTimeStepCuts="1"/>
+        maxTimeStepCuts="1"
+        maxAllowedResidualNorm="1e20"
+        lineSearchAction="None"/>
       <LinearSolverParameters
         directParallel="0"/>
     </SinglePhasePoromechanicsEmbeddedFractures>

src/coreComponents/constitutive/contact/CoulombFriction.hpp

@@ -221,7 +221,7 @@ GEOS_HOST_DEVICE
 real64 CoulombFrictionUpdates::computeLimitTangentialTractionNorm( real64 const & normalTraction,
                                                                    real64 & dLimitTangentialTractionNorm_dTraction ) const
 {
-  dLimitTangentialTractionNorm_dTraction = m_frictionCoefficient;
+  dLimitTangentialTractionNorm_dTraction = -m_frictionCoefficient;
   return ( m_cohesion - normalTraction * m_frictionCoefficient );
 }
 
@@ -274,10 +274,10 @@ inline void CoulombFrictionUpdates::computeShearTraction( localIndex const k,
 
       dTractionVector_dJump[1][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[0] / slipNorm;
       dTractionVector_dJump[1][1] = limitTau * pow( slip[1], 2 )  / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
-      dTractionVector_dJump[1][2] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
+      dTractionVector_dJump[1][2] = -limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
 
       dTractionVector_dJump[2][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[1] / slipNorm;
-      dTractionVector_dJump[2][1] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
+      dTractionVector_dJump[2][1] = -limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
       dTractionVector_dJump[2][2] = limitTau * pow( slip[0], 2 )  / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 );
 
       // Compute elastic component of the slip for this case

src/coreComponents/constitutive/contact/FrictionBase.hpp

@@ -175,7 +175,11 @@ class FrictionBaseUpdates
   inline
   virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction,
                                                      real64 & dLimitTangentialTractionNorm_dTraction ) const
-  { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0; };
+  {
+    GEOS_UNUSED_VAR( normalTraction );
+    dLimitTangentialTractionNorm_dTraction = 0.0;
+    return 0;
+  }
 
 protected:
 

src/coreComponents/constitutive/contact/FrictionlessContact.hpp

@@ -63,20 +63,6 @@ class FrictionlessContactUpdates : public FrictionBaseUpdates
                                     arraySlice1d< real64 const > const & tractionVector,
                                     integer & fractureState ) const override final;
 
-
-  /**
-   * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction
-   * @param[in] normalTraction the normal traction
-   * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction
-   * @return the limit tangential traction norm
-   */
-  GEOS_HOST_DEVICE
-  inline
-  virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction,
-                                                     real64 & dLimitTangentialTractionNorm_dTraction ) const override final
-  { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0.0; }
-
-private:
 };
 
 

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp

@@ -1593,7 +1593,7 @@ void SolidMechanicsLagrangeContact::
                   {
                     elemRHS[i] = Ja * ( traction[kfe][i] - limitTau * sliding[ i-1 ] / slidingNorm );
 
-                    dRdT( i, 0 ) = Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm;
+                    dRdT( i, 0 ) = -Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm;
                     dRdT( i, i ) = Ja;
                   }
 
@@ -1641,7 +1641,7 @@ void SolidMechanicsLagrangeContact::
                     {
                       elemRHS[i] = Ja * traction[kfe][i] * ( 1.0 - limitTau / vauxNorm );
 
-                      dRdT( i, 0 ) = Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm;
+                      dRdT( i, 0 ) = -Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm;
                       dRdT( i, i ) = Ja;
                     }
                   }

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp

@@ -471,7 +471,7 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups()
                                                                      [&]( localIndex const,
                                                                           SurfaceElementRegion & region )
     {
-      region.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion )
+      region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion )
       {
         subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ).
           setApplyDefaultValue( region.getDefaultAperture() );
@@ -691,8 +691,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_
     {
       saveConvergedState( subRegion );
 
-      arrayView1d< real64 > const & dVol = subRegion.template getField< fields::flow::deltaVolume >();
-      dVol.zero();
+      applyDeltaVolume( subRegion );
 
       // This should fix NaN density in newly created fracture elements
       updatePorosityAndPermeability( subRegion );
@@ -707,8 +706,8 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_
 
     } );
 
-    mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const,
-                                                                                          FaceElementSubRegion & subRegion )
+    mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+                                                                                             SurfaceElementSubRegion & subRegion )
     {
       arrayView1d< real64 const > const aper = subRegion.getField< fields::flow::hydraulicAperture >();
       arrayView1d< real64 > const aper0 = subRegion.getField< fields::flow::aperture0 >();
@@ -728,6 +727,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_
       fluid.saveConvergedState();
 
     } );
+
   } );
 }
 
@@ -755,14 +755,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
       singlePhaseBaseKernels::StatisticsKernel::
         saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres );
 
-      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;
-      } );
+      applyDeltaVolume( subRegion );
 
       SingleFluidBase const & fluid =
         getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
@@ -791,8 +784,8 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
 
     } );
 
-    mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const,
-                                                                                          FaceElementSubRegion & subRegion )
+    mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+                                                                                             SurfaceElementSubRegion & subRegion )
     {
       arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank();
       arrayView1d< real64 const > const volume = subRegion.getElementVolume();
@@ -817,7 +810,6 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
         }
       } );
     } );
-
   } );
 }
 
@@ -1403,4 +1395,15 @@ void SinglePhaseBase::saveConvergedState( ElementSubRegionBase & subRegion ) con
   mass_n.setValues< parallelDevicePolicy<> >( mass );
 }
 
+void SinglePhaseBase::applyDeltaVolume( ElementSubRegionBase & subRegion )
+{
+  arrayView1d< real64 > const dVol = subRegion.template getField< fields::flow::deltaVolume >();
+  arrayView1d< real64 > const vol = subRegion.template getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString());
+  forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    vol[ei] += dVol[ei];
+    dVol[ei] = 0.0;
+  } );
+}
+
 } /* namespace geos */

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp

@@ -387,6 +387,8 @@ class SinglePhaseBase : public FlowSolverBase
    */
   virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override;
 
+  void applyDeltaVolume( ElementSubRegionBase & subRegion );
+
   /**
    * @brief Structure holding views into fluid properties used by the base solver.
    */

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp

@@ -26,6 +26,7 @@
 #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp"
+#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
 #include "finiteVolume/FluxApproximationBase.hpp"
@@ -510,10 +511,10 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d
           using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel );
           typename HydraulicApertureModelType::KernelWrapper hydraulicApertureModelWrapper = castedHydraulicApertureModel.createKernelWrapper();
 
-          poromechanicsEFEMKernels::StateUpdateKernel::
+          poromechanicsFracturesKernels::StateUpdateKernel::
             launch< parallelDevicePolicy<> >( subRegion.size(),
-                                              hydraulicApertureModelWrapper,
                                               porousMaterialWrapper,
+                                              hydraulicApertureModelWrapper,
                                               dispJump,
                                               pressure,
                                               area,

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp

@@ -250,7 +250,6 @@ class SinglePhasePoromechanicsEFEM :
   /// The rank global densities
   arrayView2d< real64 const > const m_solidDensity;
   arrayView2d< real64 const > const m_fluidDensity;
-  arrayView2d< real64 const > const m_fluidDensity_n;
   arrayView2d< real64 const > const m_dFluidDensity_dPressure;
 
   /// The rank-global fluid pressure array.
@@ -259,9 +258,6 @@ class SinglePhasePoromechanicsEFEM :
   /// The rank-global fluid pressure array.
   arrayView1d< real64 const > const m_fracturePressure;
 
-  /// The rank-global delta-fluid pressure array.
-  arrayView2d< real64 const > const m_porosity_n;
-
   arrayView2d< real64 const > const m_tractionVec;
 
   arrayView3d< real64 const > const m_dTraction_dJump;
@@ -278,10 +274,14 @@ class SinglePhasePoromechanicsEFEM :
 
   arrayView1d< real64 const > const m_surfaceArea;
 
-  arrayView1d< real64 const > const m_elementVolume;
+  arrayView1d< real64 const > const m_elementVolumeCell;
+
+  arrayView1d< real64 const > const m_elementVolumeFrac;
 
   arrayView1d< real64 const > const m_deltaVolume;
 
+  arrayView1d< real64 const > const m_mass_n;
+
   SortedArrayView< localIndex const > const m_fracturedElems;
 
   ArrayOfArraysView< localIndex const > const m_cellsToEmbeddedSurfaces;
@@ -305,70 +305,6 @@ using SinglePhaseKernelFactory = finiteElement::KernelFactory< SinglePhasePorome
                                                                real64 const (&)[3],
                                                                string const >;
 
-/**
- * @brief A struct to perform volume, aperture and fracture traction updates
- */
-struct StateUpdateKernel
-{
-
-  /**
-   * @brief Launch the kernel function doing volume, aperture and fracture traction updates
-   * @tparam POLICY the type of policy used in the kernel launch
-   * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates
-   * @param[in] size the size of the subregion
-   * @param[in] contactWrapper the wrapper implementing the contact relationship
-   * @param[in] dispJump the displacement jump
-   * @param[in] pressure the pressure
-   * @param[in] area the area
-   * @param[in] volume the volume
-   * @param[out] deltaVolume the change in volume
-   * @param[out] aperture the aperture
-   * @param[out] hydraulicAperture the effecture aperture
-   * @param[out] fractureContactTraction the fracture contact traction
-   */
-  template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER >
-  static void
-  launch( localIndex const size,
-          CONTACT_WRAPPER const & contactWrapper,
-          POROUS_WRAPPER const & porousMaterialWrapper,
-          arrayView2d< real64 const > const & dispJump,
-          arrayView1d< real64 const > const & pressure,
-          arrayView1d< real64 const > const & area,
-          arrayView1d< real64 const > const & volume,
-          arrayView1d< real64 > const & deltaVolume,
-          arrayView1d< real64 > const & aperture,
-          arrayView1d< real64 const > const & oldHydraulicAperture,
-          arrayView1d< real64 > const & hydraulicAperture,
-          arrayView2d< real64 > const & fractureEffectiveTraction )
-  {
-    forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
-    {
-      // update aperture to be equal to the normal displacement jump
-      aperture[k] = dispJump[k][0]; // the first component of the jump is the normal one.
-
-      real64 dHydraulicAperture_dNormalJump = 0.0;
-      real64 dHydraulicAperture_dNormalTraction = 0.0;
-      hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k],
-                                                                      fractureEffectiveTraction[k][0],
-                                                                      dHydraulicAperture_dNormalJump,
-                                                                      dHydraulicAperture_dNormalTraction );
-
-      deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k];
-
-      real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] );
-      real64 const effectiveTraction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureEffectiveTraction[k] );
-
-      // all perm update models below should need effective traction instead of total traction
-      // (total traction is combined forces of fluid pressure and effective traction)
-      porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( k, 0, pressure[k],
-                                                                            oldHydraulicAperture[k], hydraulicAperture[k],
-                                                                            dHydraulicAperture_dNormalJump,
-                                                                            jump, effectiveTraction );
-
-    } );
-  }
-};
-
 } // namespace poromechanicsEFEMKernels
 
 } /* namespace geos */

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp

@@ -77,12 +77,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
   m_wDofNumber( jumpDofNumber ),
   m_solidDensity( inputConstitutiveType.getDensity() ),
   m_fluidDensity( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ),
-  m_fluidDensity_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ),
   m_dFluidDensity_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                      fluidModelKey ) ).dDensity_dPressure() ),
   m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ),
   m_fracturePressure( embeddedSurfSubRegion.template getField< fields::flow::pressure >() ),
-  m_porosity_n( inputConstitutiveType.getPorosity_n() ),
   m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >() ),
   m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >() ),
   m_dTraction_dPressure( embeddedSurfSubRegion.getField< fields::contact::dTraction_dPressure >() ),
@@ -91,8 +89,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
   m_tVec2( embeddedSurfSubRegion.getTangentVector2() ),
   m_surfaceCenter( embeddedSurfSubRegion.getElementCenter() ),
   m_surfaceArea( embeddedSurfSubRegion.getElementArea() ),
-  m_elementVolume( elementSubRegion.getElementVolume() ),
-  m_deltaVolume( elementSubRegion.template getField< fields::flow::deltaVolume >() ),
+  m_elementVolumeCell( elementSubRegion.getElementVolume() ),
+  m_elementVolumeFrac( embeddedSurfSubRegion.getElementVolume() ),
+  m_deltaVolume( embeddedSurfSubRegion.template getField< fields::flow::deltaVolume >() ),
+  m_mass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ),
   m_fracturedElems( elementSubRegion.fracturedElementsList() ),
   m_cellsToEmbeddedSurfaces( elementSubRegion.embeddedSurfacesList().toViewConst() ),
   m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] },
@@ -148,7 +148,7 @@ setup( localIndex const k,
 {
   localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];
 
-  stack.hInv = m_surfaceArea[embSurfIndex] / m_elementVolume[k];
+  stack.hInv = m_surfaceArea[embSurfIndex] / m_elementVolumeCell[k];
   for( localIndex a=0; a<numNodesPerElem; ++a )
   {
     localIndex const localNodeIndex = m_elemsToNodes( k, a );
@@ -320,10 +320,8 @@ complete( localIndex const k,
   real64 const localJumpFracPressureJacobian = m_surfaceArea[embSurfIndex];
 
   // Mass balance accumulation
-  real64 const newVolume = m_elementVolume( embSurfIndex ) + m_deltaVolume( embSurfIndex );
-  real64 const newMass =  m_fluidDensity( embSurfIndex, 0 ) * newVolume;
-  real64 const oldMass =  m_fluidDensity_n( embSurfIndex, 0 ) * m_elementVolume( embSurfIndex );
-  real64 const localFlowResidual = ( newMass - oldMass );
+  real64 const newVolume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex );
+  real64 const localFlowResidual = m_fluidDensity( embSurfIndex, 0 ) * newVolume - m_mass_n[embSurfIndex];
   real64 const localFlowJumpJacobian = m_fluidDensity( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ];
   real64 const localFlowFlowJacobian = m_dFluidDensity_dPressure( embSurfIndex, 0 ) * newVolume;