feat: Add reference thermal conductivity

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ 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 #3959 (2026-02-03) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3959-15395-0c3212e.tar.gz>
+=====================
+Addition of reference thermal conductivity for thermal single phase
+
 PR #3940 (2026-01-27) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3940-15307-53de7ba.tar.gz>
 =====================
 Fix the transimissibility calculated between a cell and a surface element

inputFiles/singlePhaseWell/.gitattributes

@@ -0,0 +1 @@
+mesh_1d.vtr filter=lfs diff=lfs merge=lfs -text

inputFiles/singlePhaseWell/singlePhaseWell.ats

@@ -57,6 +57,13 @@ decks = [
         partitions=((1, 1, 1), (2, 2, 2)),
         restart_step=5,
         check_step=10,
+        restartcheck_params=RestartcheckParameters(**restartcheck_params)),
+    TestDeck(
+        name="thermal_compressible_single_phase_wells_1d",
+        description='Thermal well test (1D, thermal, 1 wells)',
+        partitions=((1, 1, 1),),
+        restart_step=2,
+        check_step=0,
         restartcheck_params=RestartcheckParameters(**restartcheck_params))
 ]
 

inputFiles/singlePhaseWell/thermal_compressible_single_phase_wells_1d.xml

@@ -0,0 +1,206 @@
+<?xml version="1.0"?>
+<Problem>
+  <Solvers>
+    <SinglePhaseReservoir
+      name="RESERVOIR.SOLVER"
+      flowSolverName="FLOW.SOLVER"
+      wellSolverName="WELL.SOLVER"
+      logLevel="1"
+      initialDt="86400"
+      targetRegions="{ RESERVOIR, WELL }">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-4"
+        maxAllowedResidualNorm="1e12"
+        newtonMaxIter="25" />
+      <LinearSolverParameters
+        solverType="fgmres"
+        preconditionerType="mgr"
+        krylovTol="1e-4"
+        krylovAdaptiveTol="1"
+        krylovWeakestTol="1e-4"
+        logLevel="1" />
+    </SinglePhaseReservoir>
+
+    <SinglePhaseFVM
+      name="FLOW.SOLVER"
+      discretization="FLUID.TPFA"
+      targetRegions="{ RESERVOIR }"
+      temperature="380.15"
+      isThermal="1"
+      logLevel="2" />
+    <SinglePhaseWell
+      name="WELL.SOLVER"
+      logLevel="1"
+      isThermal="1"
+      targetRegions="{ WELL }">
+      <WellControls
+        name="WELL.CONTROL"
+        type="injector"
+        control="totalVolRate"
+        referenceElevation="-200"
+        useSurfaceConditions="1"
+        surfacePressure="101325.0"
+        surfaceTemperature="283.15"
+        injectionTemperature="293.15"
+        initialPressureCoefficient="0"
+        injectionStream="{1.0}"
+        targetBHP="1e10"
+        enableCrossflow="0"
+        targetTotalRateTableName="RATE" />
+    </SinglePhaseWell>
+  </Solvers>
+
+  <NumericalMethods>
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="FLUID.TPFA" />
+    </FiniteVolume>
+  </NumericalMethods>
+
+  <Mesh>
+    <VTKMesh
+      name="GLOBAL"
+      file="mesh_1d.vtr"
+      fieldNamesInGEOS="{ THERMAL.COND_referenceThermalConductivity }"
+      fieldsToImport="{ CONDUCTIVITY }"
+      logLevel="5">
+      <InternalWell
+        name="well.INJECTOR"
+        wellRegionName="WELL"
+        wellControlsName="WELL.CONTROL"
+        polylineNodeCoords="{
+          { 6440.37, 0.0, -200 },
+          { 6440.37, 0.0, -300 }
+        }"
+        polylineSegmentConn="{ { 0, 1 } }"
+        radius="0.19"
+        numElementsPerSegment="1">
+        <Perforation
+          name="WELL.PERF.1"
+          distanceFromHead="75.0" />
+      </InternalWell>
+    </VTKMesh>
+  </Mesh>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="RESERVOIR"
+      cellBlocks="{ hexahedra }"
+      materialList="{ FLUID, ROCK, THERMAL.COND  }" />
+    <WellElementRegion
+      name="WELL" materialList="{ FLUID }" />
+  </ElementRegions>
+
+  <Constitutive>
+    <ThermalCompressibleSinglePhaseFluid
+      name="FLUID"
+      defaultDensity="1032.26544679"
+      defaultViscosity="0.0002634"
+      referencePressure="40810068.4183"
+      referenceDensity="1032.26544679"
+      compressibility="0.00037130e-5"
+      referenceViscosity="0.0002634"
+      referenceTemperature="380.15"
+      thermalExpansionCoeff="3e-4"
+      viscosibility="0.0"
+      specificHeatCapacity="4180.0"
+      referenceInternalEnergy="0.001" />
+
+    <NullModel
+      name="NULL.SOLID" />
+
+    <SinglePhaseThermalConductivity
+      name="THERMAL.COND"
+      defaultThermalConductivityComponents="{ 0.68337, 0.68337, 0.68337 }"
+      thermalConductivityGradientComponents="{ -1.2e-5, -1.2e-5, -1.5e-5 }"
+      referenceTemperature="380.15" />
+
+    <SolidInternalEnergy
+      name="ROCK.ENERGY"
+      referenceVolumetricHeatCapacity="2.07e6"
+      referenceTemperature="273.0"
+      referenceInternalEnergy="0.0" />
+
+    <CompressibleSolidConstantPermeability
+      name="ROCK"
+      solidModelName="NULL.SOLID"
+      porosityModelName="ROCK.POROSITY"
+      permeabilityModelName="ROCK.PERMEABILITY"
+      solidInternalEnergyModelName="ROCK.ENERGY" />
+
+    <PressurePorosity
+      name="ROCK.POROSITY"
+      defaultReferencePorosity="0.3"
+      referencePressure="40800000.0000"
+      compressibility="4.0950e-14" />
+
+    <ConstantPermeability
+      name="ROCK.PERMEABILITY"
+      permeabilityComponents="{ 9.869233e-14, 9.869233e-14, 9.869233e-14}" />
+
+  </Constitutive>
+
+  <FieldSpecifications>
+    <FieldSpecification
+      name="INITIAL.PRESSURE"
+      setNames="{ all }"
+      objectPath="ElementRegions/RESERVOIR"
+      fieldName="pressure"
+      scale="169.12e5"
+      initialCondition="1" />
+    <FieldSpecification
+      name="INITIAL.TEMPERATURE"
+      setNames="{ all }"
+      objectPath="ElementRegions/RESERVOIR"
+      fieldName="temperature"
+      scale="380.15"
+      initialCondition="1" />
+  </FieldSpecifications>
+
+  <Functions>
+    <TableFunction
+      name="RATE"
+      interpolation="lower"
+      inputVarNames="{ time }"
+      coordinates="{ 0.0, 2592000, 86400000}"
+      values="{ 0.0, 2.911949431, 2.911949431 }" />
+
+  </Functions>
+
+  <Tasks>
+    <SinglePhaseStatistics
+      name="STATISTICS"
+      flowSolverName="FLOW.SOLVER"
+      logLevel="1"
+      writeCSV="1" />
+  </Tasks>
+
+  <Outputs>
+    <VTK
+      name="VTK" />
+    <Restart
+      name="RESTART"/>      
+  </Outputs>
+
+  <Events maxTime="9504000">
+    <PeriodicEvent
+      name="SOLVER.APPLICATION"
+      target="/Solvers/RESERVOIR.SOLVER"
+      maxEventDt="864000" />
+    <PeriodicEvent
+      name="STATISTICS"
+      target="/Tasks/STATISTICS"
+      targetExactTimestep="1"
+      timeFrequency="864000" />
+    <PeriodicEvent
+      name="VTK"
+      target="/Outputs/VTK"
+      targetExactTimestep="1"
+      timeFrequency="864000" />
+    <PeriodicEvent
+      name="RESTART"
+      timeFrequency="1728000"
+      targetExactTimestep="1"
+      target="/Outputs/RESTART" />
+  </Events>
+</Problem>

src/coreComponents/constitutive/thermalConductivity/SinglePhaseThermalConductivity.cpp

@@ -46,21 +46,35 @@ SinglePhaseThermalConductivity::SinglePhaseThermalConductivity( string const & n
     setDescription( "The reference temperature at which the conductivity components are equal to the default values" );
 }
 
+void SinglePhaseThermalConductivity::allocateConstitutiveData( Group & parent, localIndex const numPts )
+{
+  SinglePhaseThermalConductivityBase::allocateConstitutiveData( parent, numPts );
+
+  initializeState( parent.size());
+}
+
+void SinglePhaseThermalConductivity::initializeState() const
+{
+  SinglePhaseThermalConductivityBase::initializeState();
+  localIndex const size = m_referenceThermalConductivity.size( 0 );
+  initializeState( size );
+}
+
 void SinglePhaseThermalConductivity::initializeRockFluidState( arrayView2d< real64 const > const & initialPorosity ) const
 {
+  arrayView3d< real64 const > referenceThermalConductivity = m_referenceThermalConductivity.toViewConst();
   arrayView3d< real64 > dEffectiveConductivity_dT = m_dEffectiveConductivity_dT.toView();
   arrayView3d< real64 > effectiveConductivity = m_effectiveConductivity.toView();
-  R1Tensor const defaultThermalConductivityComponents = m_defaultThermalConductivityComponents;
   R1Tensor const thermalConductivityGradientComponents = m_thermalConductivityGradientComponents;
 
   forAll< parallelDevicePolicy<> >( initialPorosity.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const ei )
   {
     // NOTE: enforcing 1 quadrature point
     for( localIndex q = 0; q < 1; ++q )
     {
-      effectiveConductivity[ei][q][0] = defaultThermalConductivityComponents[0];
-      effectiveConductivity[ei][q][1] = defaultThermalConductivityComponents[1];
-      effectiveConductivity[ei][q][2] = defaultThermalConductivityComponents[2];
+      effectiveConductivity[ei][q][0] = referenceThermalConductivity[ei][q][0];
+      effectiveConductivity[ei][q][1] = referenceThermalConductivity[ei][q][1];
+      effectiveConductivity[ei][q][2] = referenceThermalConductivity[ei][q][2];
 
       dEffectiveConductivity_dT[ei][q][0] = thermalConductivityGradientComponents[0];
       dEffectiveConductivity_dT[ei][q][1] = thermalConductivityGradientComponents[1];
@@ -71,9 +85,9 @@ void SinglePhaseThermalConductivity::initializeRockFluidState( arrayView2d< real
 
 void SinglePhaseThermalConductivity::updateFromTemperature( arrayView1d< real64 const > const & temperature ) const
 {
+  arrayView3d< real64 const > referenceThermalConductivity = m_referenceThermalConductivity.toViewConst();
   arrayView3d< real64 > dEffectiveConductivity_dT = m_dEffectiveConductivity_dT.toView();
   arrayView3d< real64 > effectiveConductivity = m_effectiveConductivity.toView();
-  R1Tensor const defaultThermalConductivityComponents = m_defaultThermalConductivityComponents;
   R1Tensor const thermalConductivityGradientComponents = m_thermalConductivityGradientComponents;
   real64 const referenceTemperature = m_referenceTemperature;
 
@@ -84,9 +98,9 @@ void SinglePhaseThermalConductivity::updateFromTemperature( arrayView1d< real64
 
       real64 const deltaTemperature = temperature[ei] - referenceTemperature;
 
-      effectiveConductivity[ei][q][0] = defaultThermalConductivityComponents[0] + thermalConductivityGradientComponents[0] * deltaTemperature;
-      effectiveConductivity[ei][q][1] = defaultThermalConductivityComponents[1] + thermalConductivityGradientComponents[1] * deltaTemperature;
-      effectiveConductivity[ei][q][2] = defaultThermalConductivityComponents[2] + thermalConductivityGradientComponents[2] * deltaTemperature;
+      effectiveConductivity[ei][q][0] = referenceThermalConductivity[ei][q][0] + thermalConductivityGradientComponents[0] * deltaTemperature;
+      effectiveConductivity[ei][q][1] = referenceThermalConductivity[ei][q][1] + thermalConductivityGradientComponents[1] * deltaTemperature;
+      effectiveConductivity[ei][q][2] = referenceThermalConductivity[ei][q][2] + thermalConductivityGradientComponents[2] * deltaTemperature;
 
       for( localIndex i=0; i<=2; i++ )
       {
@@ -114,6 +128,24 @@ void SinglePhaseThermalConductivity::postInputInitialization()
 
 }
 
+void SinglePhaseThermalConductivity::initializeState( localIndex const size ) const
+{
+  integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point
+
+  arrayView3d< real64 > referenceThermalConductivity = m_referenceThermalConductivity.toView();
+  R1Tensor const defaultThermalConductivityComponents = m_defaultThermalConductivityComponents;
+
+  forAll< parallelDevicePolicy<> >( size, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    for( localIndex q = 0; q < numQuad; ++q )
+    {
+      referenceThermalConductivity[ei][q][0] = defaultThermalConductivityComponents[0];
+      referenceThermalConductivity[ei][q][1] = defaultThermalConductivityComponents[1];
+      referenceThermalConductivity[ei][q][2] = defaultThermalConductivityComponents[2];
+    }
+  } );
+}
+
 REGISTER_CATALOG_ENTRY( ConstitutiveBase, SinglePhaseThermalConductivity, string const &, Group * const )
 
 } // namespace constitutive

src/coreComponents/constitutive/thermalConductivity/SinglePhaseThermalConductivity.hpp

@@ -40,9 +40,11 @@ class SinglePhaseThermalConductivityUpdate : public SinglePhaseThermalConductivi
    * @param effectiveConductivity the array of cell-wise effective conductivities in the subregion
    * the subregion
    */
-  SinglePhaseThermalConductivityUpdate( arrayView3d< real64 > const & effectiveConductivity,
+  SinglePhaseThermalConductivityUpdate( arrayView3d< real64 > const & referenceThermalConductivity,
+                                        arrayView3d< real64 > const & effectiveConductivity,
                                         arrayView3d< real64 > const & dEffectiveConductivity_dT )
-    : SinglePhaseThermalConductivityBaseUpdate( effectiveConductivity,
+    : SinglePhaseThermalConductivityBaseUpdate( referenceThermalConductivity,
+                                                effectiveConductivity,
                                                 dEffectiveConductivity_dT )
   {}
 
@@ -72,10 +74,16 @@ class SinglePhaseThermalConductivity : public SinglePhaseThermalConductivityBase
 
   virtual string getCatalogName() const override { return catalogName(); }
 
+  virtual void allocateConstitutiveData( Group & parent, localIndex const numPts ) override;
+
   virtual void initializeRockFluidState( arrayView2d< real64 const > const & initialPorosity ) const override final;
 
   virtual void updateFromTemperature( arrayView1d< real64 const > const & temperature ) const override final;
 
+  virtual void initializeState() const override;
+
+  void initializeState( localIndex const size ) const;
+
   /// Type of kernel wrapper for in-kernel update
   using KernelWrapper = SinglePhaseThermalConductivityUpdate;
 
@@ -85,7 +93,8 @@ class SinglePhaseThermalConductivity : public SinglePhaseThermalConductivityBase
    */
   KernelWrapper createKernelWrapper() const
   {
-    return KernelWrapper( m_effectiveConductivity,
+    return KernelWrapper( m_referenceThermalConductivity,
+                          m_effectiveConductivity,
                           m_dEffectiveConductivity_dT );
   }
 

src/coreComponents/constitutive/thermalConductivity/SinglePhaseThermalConductivityBase.cpp

@@ -31,19 +31,24 @@ namespace constitutive
 SinglePhaseThermalConductivityBase::SinglePhaseThermalConductivityBase( string const & name, Group * const parent )
   : ConstitutiveBase( name, parent )
 {
+  registerField< fields::thermalconductivity::referenceThermalConductivity >( &m_referenceThermalConductivity );
   registerField< fields::thermalconductivity::effectiveConductivity >( &m_effectiveConductivity );
   registerField< fields::thermalconductivity::dEffectiveConductivity_dT >( &m_dEffectiveConductivity_dT );
 }
 
 void SinglePhaseThermalConductivityBase::allocateConstitutiveData( Group & parent, localIndex const numPts )
 {
   // NOTE: enforcing 1 quadrature point
+  m_referenceThermalConductivity.resize( 0, 1, 3 );
   m_effectiveConductivity.resize( 0, 1, 3 );
   m_dEffectiveConductivity_dT.resize( 0, 1, 3 );
 
   ConstitutiveBase::allocateConstitutiveData( parent, numPts );
 }
 
+void SinglePhaseThermalConductivityBase::initializeState() const
+{}
+
 } // namespace constitutive
 
 } // namespace geos