fix: delta volume fix plus some other stuff

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp

@@ -638,8 +638,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 );
@@ -703,14 +702,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() ) );
@@ -1350,4 +1342,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

@@ -391,6 +391,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"
@@ -520,10 +521,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;
@@ -284,6 +280,8 @@ class SinglePhasePoromechanicsEFEM :
 
   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;
@@ -307,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 >() ),
@@ -94,6 +92,7 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
   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] },
@@ -322,9 +321,7 @@ complete( localIndex const k,
 
   // Mass balance accumulation
   real64 const newVolume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex );
-  real64 const newMass = m_fluidDensity( embSurfIndex, 0 ) * newVolume;
-  real64 const oldMass = m_fluidDensity_n( embSurfIndex, 0 ) * m_elementVolumeFrac( embSurfIndex );
-  real64 const localFlowResidual = ( newMass - oldMass );
+  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;
 

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

@@ -57,9 +57,7 @@ class ThermalSinglePhasePoromechanicsEFEM :
   using Base::m_matrixPresDofNumber;
   using Base::m_wDofNumber;
   using Base::m_fluidDensity;
-  using Base::m_fluidDensity_n;
   using Base::m_dFluidDensity_dPressure;
-  using Base::m_porosity_n;
   using Base::m_surfaceArea;
   using Base::m_elementVolumeFrac;
   using Base::m_deltaVolume;
@@ -164,17 +162,18 @@ class ThermalSinglePhasePoromechanicsEFEM :
   arrayView2d< real64 const > const m_dFluidDensity_dTemperature;
 
   /// Views on fluid internal energy
-  arrayView2d< real64 const > const m_fluidInternalEnergy_n;
   arrayView2d< real64 const > const m_fluidInternalEnergy;
   arrayView2d< real64 const > const m_dFluidInternalEnergy_dPressure;
   arrayView2d< real64 const > const m_dFluidInternalEnergy_dTemperature;
 
   /// Views on temperature
-  arrayView1d< real64 const > const m_temperature_n;
   arrayView1d< real64 const > const m_temperature;
 
   /// The rank-global fluid pressure array.
   arrayView1d< real64 const > const m_matrixTemperature;
+
+  /// Views on energy
+  arrayView1d< real64 const > const m_energy_n;
 };
 
 

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

@@ -68,15 +68,14 @@ ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager,
         fluidModelKey ),
   m_dFluidDensity_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                         fluidModelKey ) ).dDensity_dTemperature() ),
-  m_fluidInternalEnergy_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ),
   m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ),
   m_dFluidInternalEnergy_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                             fluidModelKey ) ).dInternalEnergy_dPressure() ),
   m_dFluidInternalEnergy_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                                fluidModelKey ) ).dInternalEnergy_dTemperature() ),
-  m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ),
   m_temperature( embeddedSurfSubRegion.template getField< fields::flow::temperature >() ),
-  m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() )
+  m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ),
+  m_energy_n( embeddedSurfSubRegion.template getField< fields::flow::energy_n >() )
 
 {}
 
@@ -176,9 +175,8 @@ complete( localIndex const k,
   localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];
   // Energy balance accumulation
   real64 const volume        =  m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex );
-  real64 const volume_n      =  m_elementVolumeFrac( embSurfIndex );
   real64 const fluidEnergy   =  m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume;
-  real64 const fluidEnergy_n =  m_fluidDensity_n( embSurfIndex, 0 ) * m_fluidInternalEnergy_n( embSurfIndex, 0 ) * volume_n;
+  real64 const fluidEnergy_n =  m_energy_n[embSurfIndex];
 
   stack.dFluidMassIncrement_dTemperature =  m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * volume;