Surface generator fix for diffusion

src/coreComponents/constitutive/diffusion/ConstantDiffusion.cpp

@@ -52,6 +52,27 @@ void ConstantDiffusion::allocateConstitutiveData( Group & parent, localIndex con
   }
 }
 
+void ConstantDiffusion::initializeTemperatureState( arrayView1d< real64 const > const & initialTemperature ) const
+{
+  DiffusionBase::initializeTemperatureState( initialTemperature );
+
+  localIndex const numE = m_diffusivity.size( 0 );
+  integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point
+
+  auto diffusivityView = m_diffusivity.toView();
+  auto const diffusivityComponentsView = m_diffusivityComponents.toViewConst();
+
+  forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    for( localIndex q = 0; q < numQuad; ++q )
+    {
+      diffusivityView[ei][q][0] = diffusivityComponentsView[0];
+      diffusivityView[ei][q][1] = diffusivityComponentsView[1];
+      diffusivityView[ei][q][2] = diffusivityComponentsView[2];
+    }
+  } );
+}
+
 void ConstantDiffusion::postInputInitialization()
 {
   GEOS_THROW_IF( m_diffusivityComponents.size() != 3,

src/coreComponents/constitutive/diffusion/ConstantDiffusion.hpp

@@ -73,6 +73,8 @@ class ConstantDiffusion : public DiffusionBase
   virtual void allocateConstitutiveData( dataRepository::Group & parent,
                                          localIndex const numPts ) override final;
 
+  virtual void initializeTemperatureState( arrayView1d< real64 const > const & initialTemperature ) const override;
+
   /// Type of kernel wrapper for in-kernel update
   using KernelWrapper = ConstantDiffusionUpdate;
 

src/coreComponents/constitutive/diffusion/DiffusionBase.cpp

@@ -88,6 +88,29 @@ void DiffusionBase::allocateConstitutiveData( Group & parent, localIndex const n
   }
 }
 
+void DiffusionBase::initializeTemperatureState( arrayView1d< real64 const > const & initialTemperature ) const
+{
+  GEOS_UNUSED_VAR( initialTemperature );
+
+  localIndex const numE = m_phaseDiffusivityMultiplier.size( 0 );
+  integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point
+  integer const numPhases = numFluidPhases();
+
+  auto phaseDiffusivityMultiplierView = m_phaseDiffusivityMultiplier.toView();
+  auto const defaultPhaseDiffusivityMultiplierView = m_defaultPhaseDiffusivityMultiplier.toViewConst();
+
+  forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    for( localIndex q = 0; q < numQuad; ++q )
+    {
+      for( integer ip = 0; ip < numPhases; ++ip )
+      {
+        phaseDiffusivityMultiplierView[ei][q][ip] = defaultPhaseDiffusivityMultiplierView[ip];
+      }
+    }
+  } );
+}
+
 } // namespace constitutive
 
 } // namespace geos

src/coreComponents/constitutive/diffusion/DiffusionBase.hpp

@@ -140,8 +140,7 @@ class DiffusionBase : public ConstitutiveBase
    *
    * Note: this is needed because for now, the temperature is treated **explicitly** in the diffusion tensor
    */
-  virtual void initializeTemperatureState( arrayView1d< real64 const > const & initialTemperature ) const
-  { GEOS_UNUSED_VAR( initialTemperature ); }
+  virtual void initializeTemperatureState( arrayView1d< real64 const > const & initialTemperature ) const;
 
   /**
    * @brief Save the temperature state (needed when diffusion depends on temperature)

src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp

@@ -29,6 +29,7 @@
 #include "mesh/SurfaceElementRegion.hpp"
 #include "mesh/utilities/ComputationalGeometry.hpp"
 #include "constitutive/thermalConductivity/SinglePhaseThermalConductivityBase.hpp"
+#include "constitutive/diffusion/DiffusionBase.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp"
 #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsLagrangianFEMKernels.hpp"
@@ -583,6 +584,14 @@ real64 SurfaceGenerator::solverStep( real64 const & time_n,
       SinglePhaseThermalConductivityBase & thermalCondModel = getConstitutiveModel< SinglePhaseThermalConductivityBase >( fractureSubRegion, thermalCondModelName );
       thermalCondModel.initializeState();
     }
+
+    string const diffusionModelName = getConstitutiveName< DiffusionBase >( fractureSubRegion );
+    if( !diffusionModelName.empty() )
+    {
+      DiffusionBase & diffusionModel = getConstitutiveModel< DiffusionBase >( fractureSubRegion, diffusionModelName );
+      arrayView1d< real64 const > const temperature = fractureSubRegion.template getField< flow::temperature >();
+      diffusionModel.initializeTemperatureState( temperature );
+    }
   } );
 
   return rval;
@@ -4775,6 +4784,8 @@ SurfaceGenerator::calculateRuptureRate( SurfaceElementRegion & faceElementRegion
 {
   GEOS_MARK_FUNCTION;
 
+  real64 constexpr timeTolerance = 1.0e-14;
+
   real64 maxRuptureRate = 0;
   FaceElementSubRegion & subRegion = faceElementRegion.getSubRegion< FaceElementSubRegion >( 0 );
 
@@ -4796,8 +4807,8 @@ SurfaceGenerator::calculateRuptureRate( SurfaceElementRegion & faceElementRegion
         localIndex const faceElem1 = fractureConnectorEdgesToFaceElements( kfc, kfe1 );
         if( !( m_faceElemsRupturedThisSolve.count( faceElem0 ) && m_faceElemsRupturedThisSolve.count( faceElem1 ) ) )
         {
-          real64 const deltaRuptureTime = fabs( ruptureTime( faceElem0 ) - ruptureTime( faceElem1 ));
-          if( deltaRuptureTime > 1.0e-14 * (ruptureTime( faceElem0 ) + ruptureTime( faceElem1 )) )
+          real64 const deltaRuptureTime = LvArray::math::abs( ruptureTime( faceElem0 ) - ruptureTime( faceElem1 ) );
+          if( deltaRuptureTime > timeTolerance )
           {
             real64 distance[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( elemCenter[ faceElem0 ] );
             LvArray::tensorOps::subtract< 3 >( distance, elemCenter[ faceElem1 ] );