feat: Add single phase viscosity dependency on temperature

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3959-16478-faf1698
+  baseline: integratedTests/baseline_integratedTests-pr4029-16507-7e3d8b5
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -5,6 +5,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 #4029 (2026-05-06) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr4029-16507-7e3d8b5.tar.gz>
+Add single phase viscosity dependency on temperature
+
 PR #3959 (2026-05-04) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3959-16478-faf1698.tar.gz>
 Add reference thermal conductivity
 

inputFiles/singlePhaseFlow/thermalSinglePhaseFlow.ats

@@ -45,6 +45,13 @@ decks = [
         restart_step=0,
         check_step=1,
         restartcheck_params=RestartcheckParameters(**restartcheck_params)), 
+    TestDeck(
+        name="thermalViscosity_1d",
+        description='Thermal problem with temperature dependent viscosity',
+        partitions=((1, 1, 1),),
+        restart_step=0,
+        check_step=1,
+        restartcheck_params=RestartcheckParameters(**restartcheck_params)),
     TestDeck(
         name="thermalCompressible_referenceConductivity",
         description='Thermal reference conductivity (1D, thermal)',

inputFiles/singlePhaseFlow/thermalViscosity_1d.xml

@@ -0,0 +1,183 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Solvers>
+    <SinglePhaseFVM
+      name="singleFlow"
+      logLevel="1"
+      discretization="tpfaFlow"
+      temperature="0"
+      isThermal="1"
+      targetRegions="{ region }">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-6"
+        newtonMaxIter="100"/>
+      <LinearSolverParameters
+        directParallel="0"/>
+    </SinglePhaseFVM>
+  </Solvers>
+
+  <NumericalMethods>
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="tpfaFlow"/>
+    </FiniteVolume>
+  </NumericalMethods>
+
+  <Mesh>
+    <InternalMesh
+      name="mesh"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 0, 20 }"
+      yCoords="{ 0, 1 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 20 }"
+      ny="{ 1 }"
+      nz="{ 1 }"
+      cellBlockNames="{ region }">
+    </InternalMesh>
+  </Mesh>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="region"
+      cellBlocks="{ region }"
+      materialList="{ rock, fluid, thermalCond }"/>
+  </ElementRegions>
+
+  <Constitutive>
+
+    <CompressibleSolidConstantPermeability
+      name="rock"
+      solidModelName="nullSolid"
+      porosityModelName="rockPorosity"
+      permeabilityModelName="rockPerm"
+      solidInternalEnergyModelName="rockInternalEnergy" />
+
+    <NullModel
+      name="nullSolid"/>
+
+    <PressurePorosity
+      name="rockPorosity"
+      defaultReferencePorosity="0.1"
+      referencePressure="0.0"
+      compressibility="3.0e-11"/>
+
+	<!-- SPHINX_SolidInternalEnergy_linear -->
+    <SolidInternalEnergy
+      name="rockInternalEnergy"
+      referenceVolumetricHeatCapacity="1.0e6"
+      referenceTemperature="0"
+      referenceInternalEnergy="0"/>
+	<!-- SPHINX_SolidInternalEnergy_linearEnd -->
+
+    <ConstantPermeability
+      name="rockPerm"
+      permeabilityComponents="{ 1.0e-18, 1.0e-18, 1.0e-18 }"/>
+
+	<!-- SPHINX_ThermalCompressibleSinglePhaseFluid -->
+    <ThermalCompressibleSinglePhaseFluid
+      name="fluid"
+      defaultDensity="1000"
+      defaultViscosity="0.001"
+      referencePressure="0.0"
+      referenceTemperature="293.15"
+      compressibility="5e-10"
+      thermalExpansionCoeff="3e-4"
+      viscosibility="0.0"
+      temperatureViscosityCoefficient="6.21e-4"
+      specificHeatCapacity="1"
+      referenceInternalEnergy="0.99"/>
+	<!-- SPHINX_ThermalCompressibleSinglePhaseFluidEnd -->
+
+    <!-- SPHINX_SinglePhaseThermalConductivity_linear -->
+    <SinglePhaseThermalConductivity
+      name="thermalCond"
+      defaultThermalConductivityComponents="{ 1.66, 1.66, 1.66 }"
+      thermalConductivityGradientComponents="{ 0, 0, 0 }"
+      referenceTemperature="293.15" />
+    <!-- SPHINX_SinglePhaseThermalConductivity_linearEnd -->    
+  </Constitutive>
+
+  <FieldSpecifications>
+
+	<!-- SPHINX_FieldSpecificationImposedPressure -->
+    <FieldSpecification
+      name="initialPressure"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region"
+      fieldName="pressure"
+      scale="0e6"/>
+
+    <FieldSpecification
+      name="sinkPressure"
+      setNames="{ xpos }"       
+      objectPath="faceManager"
+      fieldName="pressure"
+      scale="1e6"/>
+
+    <FieldSpecification
+      name="sourcePressure"
+      setNames="{ xneg }"     
+      objectPath="faceManager"
+      fieldName="pressure"
+      scale="-1e6"/>
+	<!-- SPHINX_FieldSpecificationImposedPressureEnd -->
+
+	<!-- SPHINX_FieldSpecificationImposedTemperature -->
+    <FieldSpecification
+      name="initialTemperature"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region"
+      fieldName="temperature"
+      scale="373.15"/>
+
+    <FieldSpecification
+      name="sinkTemperature"
+      setNames="{ xpos }"       
+      objectPath="faceManager"
+      fieldName="temperature"
+      scale="373.15"/>
+
+    <FieldSpecification
+      name="sourceTemperature"
+      setNames="{ xneg }"     
+      objectPath="faceManager"
+      fieldName="temperature"
+      scale="293.15"/>
+	<!-- SPHINX_FieldSpecificationImposedTemperatureEnd -->
+
+  </FieldSpecifications>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput"/>
+
+    <Restart
+      name="restartOutput"/>
+  </Outputs>
+
+  <Events 
+    maxTime="8640000">
+
+    <PeriodicEvent
+      name="outputs"
+      cycleFrequency="5"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplications"
+      beginTime="0"
+      maxEventDt="864000"
+      target="/Solvers/singleFlow"/>
+
+    <PeriodicEvent
+      name="restarts"
+      cycleFrequency="5"
+      target="/Outputs/restartOutput"/>
+
+  </Events>
+
+</Problem>

src/coreComponents/constitutive/docs/CompressibleSinglePhaseFluid.rst

@@ -9,31 +9,46 @@ Overview
 
 This model represents a compressible single-phase fluid with constant compressibility
 and pressure-dependent viscosity.
+For thermal simulations, fluid density and viscosity are temperature-dependent, governed by a constant thermal
+expansion coefficient and a constant temperature-viscosity coefficient, respectively.
 These assumptions are valid for slightly compressible fluids, such as water, and some
 types of oil with negligible amounts of dissolved gas.
 
 Specifically, fluid density is computed as
 
 .. math::
-   \rho(p) = \rho_0 ( e^{c_\rho(p - p_0)} + e^{\beta_f (T - T_0)} )
+   \rho(p) = \rho_0  e^{c_\rho(p - p_0)}
 
-where :math:`c_\rho` is compressibility, :math:`p_0` is reference pressure, :math:`\rho_0` is
-density at reference pressure, :math:`\beta_f` is the fluid thermal expansion coefficient, :math:`T_0` is reference temperature.
+for isothermal cases and,
+
+.. math::
+   \rho(p,T) = \rho_0  e^{c_\rho(p - p_0)}  e^{-\beta_\rho (T - T_0)}
+
+for thermal cases.
+
+where :math:`c_\rho` is compressibility, :math:`p_0` is reference pressure, :math:`\rho_0` is the
+density at the reference pressure, :math:`\beta_\rho` is the fluid thermal expansion coefficient,
+:math:`T_0` is reference temperature.
 
 Similarly,
 
 .. math::
    \mu(p) = \mu_0 e^{c_\mu(p - p_0)}
 
-where :math:`c_\mu` is viscosibility (viscosity compressibility), :math:`\mu_0` is reference viscosity.
+for isothermal cases and,
 
-Either exponent may be approximated by linear (default) or quadratic terms of Taylor series expansion.
-Currently there is no temperature dependence in the model, although it may be added in future.
+.. math::
+   \mu(p,T) = \mu_0 e^{c_\mu(p - p_0)} e^{-\beta_\mu (T - T_0)}
+
+for thermal cases.
+
+where :math:`c_\mu` is viscosibility (viscosity compressibility), :math:`\mu_0` is reference viscosity,
+:math:`\beta_\rho` is the fluid viscosity thermal coefficient.
 
 Parameters
 =========================
 
-The model is represented by ``<CompressibleSinglePhaseFluid>`` node in the input.
+For the isothermal case, the model is represented by ``<CompressibleSinglePhaseFluid>`` node in the input.
 
 The following attributes are supported:
 
@@ -51,4 +66,29 @@ Example
                                   compressibility="1e-19"
                                   referenceViscosity="0.001"
                                   viscosibility="0.0"/>
-  </Constitutive>
+  </Constitutive>
+
+For the thermal case, the model is represented by ``<ThermalCompressibleSinglePhaseFluid>`` node in the input.
+
+The following attributes are supported:
+
+.. include:: /docs/sphinx/datastructure/ThermalCompressibleSinglePhaseFluid.rst
+
+Example
+=========================
+
+.. code-block:: xml
+
+  <Constitutive>
+    <ThermalCompressibleSinglePhaseFluid  name="fluid"
+                                          defaultDensity="1000"
+                                          defaultViscosity="0.001"
+                                          referencePressure="0.0"
+                                          referenceTemperature="0"
+                                          compressibility="5e-10"
+                                          thermalExpansionCoeff="3e-4"
+                                          viscosibility="0.0"
+                                          viscosityExpansivity="6.21e-4"
+                                          specificHeatCapacity="1"
+                                          referenceInternalEnergy="0.99" />
+  </Constitutive>

src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp

@@ -32,12 +32,18 @@ namespace constitutive
 ThermalCompressibleSinglePhaseFluid::ThermalCompressibleSinglePhaseFluid( string const & name, Group * const parent ):
   CompressibleSinglePhaseFluid( name, parent )
 {
-  m_numDOF=2;
+  m_numDOF = 2;
+
   registerWrapper( viewKeyStruct::thermalExpansionCoeffString(), &m_thermalExpansionCoeff ).
     setApplyDefaultValue( 0.0 ).
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( "Fluid thermal expansion coefficient. Unit: 1/K" );
 
+  registerWrapper( viewKeyStruct::temperatureViscosityCoefficient(), &m_temperatureViscosityCoefficient ).
+    setApplyDefaultValue( 0.0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Fluid viscosity thermal expansion coefficient at the reference temperature. Unit: 1/K" );
+
   registerWrapper( viewKeyStruct::specificHeatCapacityString(), &m_specificHeatCapacity ).
     setApplyDefaultValue( 0.0 ).
     setInputFlag( InputFlags::OPTIONAL ).
@@ -80,6 +86,7 @@ void ThermalCompressibleSinglePhaseFluid::postInputInitialization()
   };
 
   checkNonnegative( m_thermalExpansionCoeff, viewKeyStruct::thermalExpansionCoeffString() );
+  checkNonnegative( m_temperatureViscosityCoefficient, viewKeyStruct::temperatureViscosityCoefficient() );
   checkNonnegative( m_specificHeatCapacity, viewKeyStruct::specificHeatCapacityString() );
   checkNonnegative( m_referenceInternalEnergy, viewKeyStruct::referenceInternalEnergyString() );
 
@@ -98,7 +105,7 @@ ThermalCompressibleSinglePhaseFluid::KernelWrapper
 ThermalCompressibleSinglePhaseFluid::createKernelWrapper()
 {
   return KernelWrapper( KernelWrapper::DensRelationType( m_referencePressure, m_referenceTemperature, m_referenceDensity, m_compressibility, -m_thermalExpansionCoeff ),
-                        KernelWrapper::ViscRelationType( m_referencePressure, m_referenceViscosity, m_viscosibility ),
+                        KernelWrapper::ViscRelationType( m_referencePressure, m_referenceTemperature, m_referenceViscosity, m_viscosibility, -m_temperatureViscosityCoefficient ),
                         KernelWrapper::IntEnergyRelationType( m_referenceTemperature, m_referenceInternalEnergy, m_specificHeatCapacity/m_referenceInternalEnergy ),
                         m_density.value,
                         m_density.derivs,

src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp

@@ -43,7 +43,7 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
 public:
 
   using DensRelationType      = ExponentialRelation< real64, DENS_EAT, 3 >;
-  using ViscRelationType      = ExponentialRelation< real64, VISC_EAT >;
+  using ViscRelationType      = ExponentialRelation< real64, VISC_EAT, 3 >;
   using IntEnergyRelationType = ExponentialRelation< real64, INTENERGY_EAT >;
   using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
 
@@ -114,19 +114,17 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
                         real64 & dEnthalpy_dPressure,
                         real64 & dEnthalpy_dTemperature ) const override
   {
-    m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure );
-    dViscosity_dTemperature = 0.0;
-
     m_densRelation.compute( pressure, temperature, density, dDensity_dPressure, dDensity_dTemperature );
 
+    m_viscRelation.compute( pressure, temperature, viscosity, dViscosity_dPressure, dViscosity_dTemperature );
+
     /// Compute the internal energy (only sensitive to temperature)
     m_intEnergyRelation.compute( temperature, internalEnergy, dInternalEnergy_dTemperature );
     dInternalEnergy_dPressure = 0.0;
 
     enthalpy = internalEnergy - m_refIntEnergy;
     dEnthalpy_dPressure = 0.0;
     dEnthalpy_dTemperature = dInternalEnergy_dTemperature;
-
   }
 
   GEOS_HOST_DEVICE
@@ -240,6 +238,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid
   struct viewKeyStruct : public CompressibleSinglePhaseFluid::viewKeyStruct
   {
     static constexpr char const * thermalExpansionCoeffString() { return "thermalExpansionCoeff"; }
+    static constexpr char const * temperatureViscosityCoefficient() { return "temperatureViscosityCoefficient"; }
     static constexpr char const * specificHeatCapacityString() { return "specificHeatCapacity"; }
     static constexpr char const * referenceTemperatureString() { return "referenceTemperature"; }
     static constexpr char const * referenceInternalEnergyString() { return "referenceInternalEnergy"; }
@@ -255,6 +254,9 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid
   /// scalar fluid thermal expansion coefficient
   real64 m_thermalExpansionCoeff;
 
+  /// scalar fluid temperature-viscosity coefficient
+  real64 m_temperatureViscosityCoefficient;
+
   /// scalar fluid volumetric heat capacity coefficient
   real64 m_specificHeatCapacity;