Reference states

src/porepy/__init__.py

@@ -159,7 +159,7 @@
 from porepy.numerics import ad
 from porepy.numerics.ad.operators import wrap_as_dense_ad_array, wrap_as_sparse_ad_array
 from porepy.numerics.ad.equation_system import EquationSystem
-from porepy.numerics.ad.ad_utils import (
+from porepy.numerics.ad.get_set_values import (
     get_solution_values,
     set_solution_values,
     shift_solution_values,

src/porepy/compositional/materials.py

@@ -363,6 +363,7 @@ class SolidConstants(Constants):
             "permeability": "m^2",
             "porosity": "-",
             "residual_aperture": "m",
+            "rock_compressibility": "Pa^-1",
             "shear_modulus": "Pa",
             "skin_factor": "-",
             "specific_heat_capacity": "J * kg^-1 * K^-1",
@@ -408,6 +409,8 @@ class SolidConstants(Constants):
 
     residual_aperture: number = 0.1
 
+    rock_compressibility: number = 0.0
+
     shear_modulus: number = 1.0
 
     skin_factor: number = 0.0

src/porepy/examples/example_params.py

@@ -37,6 +37,7 @@
         "fluid": pp.FluidComponent(**water),  # type: ignore[arg-type]
         "numerical": pp.NumericalConstants(**numerical_values),  # type: ignore[arg-type]
     },
+    "initialize_operator_reference_from_initial_values": True,
     # Meshing
     "grid_type": "cartesian",
     # Depending on the grid type, some subset of the meshing arguments below are used.

src/porepy/models/abstract_equations.py

@@ -506,20 +506,30 @@ def create_variables(self) -> None:
         """
         pass
 
-    def perturbation_from_reference(self, name: str, grids: list[pp.Grid]):
-        """Perturbation of some quantity ``name`` from its reference value.
+    def perturbation_from_thermodynamic_state(self, name: str, grids: list[pp.Grid]):
+        """Perturbation of some quantity ``name`` from its thermodynamic state.
 
         The parameter ``name`` should be the name of a mixed-in method, returning an
         AD operator for given ``grids``.
 
-        ``name`` should also be defined in the model's :attr:`reference_values`.
+        ``name`` should also be defined in the model's reference thermodynamic values.
 
         This method calls the model method with given ``name`` on given ``grids`` to
-        create an operator ``A``. It then fetches the respective reference value and
-        wraps it into an AD scalar ``A_0``. The return value is an operator ``A - A_0``.
+        create an operator ``A``. It then fetches the respective thermodynamic state
+        value and wraps it into an AD scalar ``A_0``. The return value is an operator
+        ``A - A_0``.
+
+        Note:
+            This method is for scalar constitutive baselines (thermodynamic state).
+            It is distinct from
+            :meth:`porepy.numerics.ad.operators.Operator.perturbation_from_reference`,
+            which is the AD linearization reference and may be spatially heterogeneous.
+            Use the operator-level method when the reference state is stored in the
+            equation system.
 
         Parameters:
-            name: Name of the quantity to be perturbed from a reference value.
+            name: Name of the quantity to be perturbed from a thermodynamic state
+                value.
             grids: List of subdomain or interface grids on which the quantity is
                 defined.
 
@@ -530,7 +540,7 @@ def perturbation_from_reference(self, name: str, grids: list[pp.Grid]):
         quantity = getattr(self, name)
         # This will throw an error if the attribute is not callable
         quantity_op = cast(pp.ad.Operator, quantity(grids))
-        # the reference values are a data class instance storing only numbers
+        # The thermodynamic state values are stored as scalar constants.
         quantity_ref = cast(pp.number, getattr(self.reference_variable_values, name))
         # The casting reflects the expected outcome, and is used to help linters find
         # the set_name method

src/porepy/models/constitutive_laws.py

@@ -2839,7 +2839,9 @@ def solid_enthalpy(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
         """
         c = self.solid_specific_heat_capacity(subdomains)
-        enthalpy = c * self.perturbation_from_reference("temperature", subdomains)
+        enthalpy = c * self.perturbation_from_thermodynamic_state(
+            "temperature", subdomains
+        )
         enthalpy.set_name("solid_enthalpy")
         return enthalpy
 
@@ -3561,9 +3563,10 @@ def pressure_stress(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         # The stress is simply found by the scalar_gradient operator, multiplied with
         # the pressure perturbation. The reference pressure is only defined on
         # sd_primary, thus there is no need for a subdomain projection.
-        stress: pp.ad.Operator = discr.scalar_gradient(
-            self.darcy_keyword
-        ) @ self.perturbation_from_reference("pressure", subdomains)
+        stress: pp.ad.Operator = (
+            discr.scalar_gradient(self.darcy_keyword)
+            @ self.pressure(subdomains).perturbation_from_reference()
+        )
         stress.set_name("pressure_stress")
         return stress
 
@@ -3709,9 +3712,10 @@ def thermal_stress(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
                 raise ValueError("Subdomains must be of dimension nd.")
 
         discr = pp.ad.BiotAd(self.stress_keyword, subdomains)
-        stress: pp.ad.Operator = discr.scalar_gradient(
-            self.enthalpy_keyword
-        ) @ self.perturbation_from_reference("temperature", subdomains)
+        stress: pp.ad.Operator = (
+            discr.scalar_gradient(self.enthalpy_keyword)
+            @ self.temperature(subdomains).perturbation_from_reference()
+        )
         stress.set_name("thermal_stress")
         return stress
 
@@ -4915,7 +4919,7 @@ def porosity_change_from_pressure(
         bulk_modulus = self.bulk_modulus(subdomains)
 
         # Pressure changes
-        dp = self.perturbation_from_reference("pressure", subdomains)
+        dp = self.pressure(subdomains).perturbation_from_reference()
 
         # Compute 1/N as defined in Coussy, 2004, https://doi.org/10.1002/0470092718.
         n_inv = (alpha - phi_ref) * (Scalar(1) - alpha) / bulk_modulus
@@ -5036,7 +5040,10 @@ def _mpsa_consistency(
         # The consistency is based on perturbation. If the variable is used directly,
         # results will not match if the reference state is not zero, see
         # :func:`test_without_fracture` in test_poromechanics.py.
-        dp = self.perturbation_from_reference(variable_name, subdomains)
+        variable = cast(
+            Callable[[list[pp.Grid]], pp.ad.Operator], getattr(self, variable_name)
+        )
+        dp = variable(subdomains).perturbation_from_reference()
         consistency_integrated = discr.consistency(physics_name) @ dp
 
         # Divide by cell volumes to counteract integration.
@@ -5066,6 +5073,11 @@ class BiotPoroMechanicsPorosity(pp.PorePyModel):
     """Specific storage. Normally defined in a mixin instance of
     :class:`~porepy.models.constitutive_laws.SpecificStorage`.
 
+    """
+    pressure: Callable[[pp.SubdomainsOrBoundaries], pp.ad.Operator]
+    """Pressure variable. Normally defined in a mixin instance of
+    :class:`~porepy.models.fluid_mass_balance.VariablesSinglePhaseFlow`.
+
     """
 
     def porosity_change_from_pressure(
@@ -5081,7 +5093,7 @@ def porosity_change_from_pressure(
 
         """
         specific_storage = self.specific_storage(subdomains)
-        dp = self.perturbation_from_reference("pressure", subdomains)
+        dp = self.pressure(subdomains).perturbation_from_reference()
 
         # Pressure change contribution
         pressure_contribution = specific_storage * dp
@@ -5107,9 +5119,14 @@ class ThermoPoroMechanicsPorosity(PoroMechanicsPorosity):
     """Biot coefficient. Normally defined in a mixin instance of
     :class:`~porepy.models.constitutive_laws.BiotCoefficient`.
 
+    """
+    temperature: Callable[[pp.SubdomainsOrBoundaries], pp.ad.Operator]
+    """Temperature variable. Normally defined in a mixin instance of
+    :class:`~porepy.models.energy_balance.VariablesEnergyBalance`.
+
     """
     temperature_variable: str
-    """Name of the pressure variable. Normally set by a mixin instance of
+    """Name of the temperature variable. Normally set by a mixin instance of
     :class:`~porepy.models.energy_balance.SolutionStrategyEnergyBalance`.
     """
 
@@ -5147,7 +5164,7 @@ def porosity_change_from_temperature(
         """
         if not all([sd.dim == self.nd for sd in subdomains]):
             raise ValueError("Subdomains must be of dimension nd.")
-        dtemperature = self.perturbation_from_reference("temperature", subdomains)
+        dtemperature = self.temperature(subdomains).perturbation_from_reference()
         phi_ref = self.reference_porosity(subdomains)
         beta = self.solid_thermal_expansion_coefficient(subdomains)
         alpha = self.biot_coefficient(subdomains)

src/porepy/models/fluid_property_library.py

@@ -107,6 +107,8 @@ def pressure_exponential(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         Extracted as a separate method to allow for easier combination with temperature
         dependent fluid density.
 
+        The perturbation is measured relative to the thermodynamic state values.
+
         Parameters:
             subdomains: List of subdomain grids.
 
@@ -118,7 +120,7 @@ def pressure_exponential(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
         # Reference variables are defined in a variables class which is assumed to be
         # available by mixin.
-        dp = self.perturbation_from_reference("pressure", subdomains)
+        dp = self.perturbation_from_thermodynamic_state("pressure", subdomains)
 
         # Wrap compressibility from fluid class as matrix (left multiplication with dp).
         c = self.fluid_compressibility(subdomains)
@@ -175,6 +177,8 @@ def temperature_exponential(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         Extracted as a separate method to allow for easier combination with temperature
         dependent fluid density.
 
+        The perturbation is measured relative to the thermodynamic state values.
+
         Parameters:
             subdomains: List of subdomain grids.
 
@@ -186,7 +190,7 @@ def temperature_exponential(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
         # Reference variables are defined in a variables class which is assumed to be
         # available by mixin.
-        dtemp = self.perturbation_from_reference("temperature", subdomains)
+        dtemp = self.perturbation_from_thermodynamic_state("temperature", subdomains)
         c = self.fluid_thermal_expansion(subdomains)
         return exp(Scalar(-1) * c * dtemp)
 
@@ -1369,7 +1373,7 @@ def specific_enthalpy_of_phase(self, phase: pp.Phase) -> ExtendedDomainFunctionT
 
         def h(domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
             c = self.fluid_specific_heat_capacity(cast(list[pp.Grid], domains))
-            enthalpy = c * self.perturbation_from_reference(
+            enthalpy = c * self.perturbation_from_thermodynamic_state(
                 "temperature", cast(list[pp.Grid], domains)
             )
             enthalpy.set_name("fluid_enthalpy")

src/porepy/models/protocol.py

@@ -809,23 +809,28 @@ class VariableProtocol(Protocol):
         """This protocol provides the declarations of the methods and the properties,
         typically defined in VariableMixin."""
 
-        def perturbation_from_reference(
+        def perturbation_from_thermodynamic_state(
             self, variable_name: str, grids: list[pp.Grid]
         ) -> pp.ad.Operator:
-            """Perturbation of some quantity ``name`` from its reference value.
+            """Perturbation of some quantity ``name`` from its thermodynamic state.
 
             The parameter ``name`` should be the name of a mixed-in method, returning an
             AD operator for given ``grids``.
 
-            ``name`` should also be defined in the model's :attr:`reference_values`.
+            ``name`` should also be defined in the model's reference thermodynamic
+            values.
 
             This method calls the model method with given ``name`` on given ``grids`` to
-            create an operator ``A``. It then fetches the respective reference value and
-            wraps it into an AD scalar ``A_0``. The return value is an operator
-            ``A - A_0``.
+            create an operator ``A``. It then fetches the respective thermodynamic
+            state value and wraps it into an AD scalar ``A_0``. The return value is an
+            operator ``A - A_0``.
+
+            Use the operator-level perturbation method when the reference state is
+            stored in the equation system.
 
             Parameters:
-                name: Name of the quantity to be perturbed from a reference value.
+                name: Name of the quantity to be perturbed from a thermodynamic state
+                    value.
                 grids: List of subdomain or interface grids on which the quantity is
                     defined.
 

src/porepy/models/solution_strategy.py

@@ -160,6 +160,7 @@ def prepare_simulation(self) -> None:
         # opposed to e.g. pressure or temperature.
         self.assign_thermodynamic_properties_to_phases()
         self.initial_condition()
+        self.initialize_operator_reference_values_from_initial_state()
         self.initialize_previous_iterate_and_time_step_values()
 
         # Initialize time dependent ad arrays, including those for boundary values.
@@ -196,6 +197,49 @@ def initialize_previous_iterate_and_time_step_values(self) -> None:
                 time_step_index=time_step_index,
             )
 
+    def initialize_operator_reference_values_from_initial_state(self) -> None:
+        """Initialize AD operator reference values from iterate-0 state.
+
+        This compatibility step aligns linearization references used by
+        :meth:`porepy.numerics.ad.operators.Operator.perturbation_from_reference`
+        with initialized primary-variable values for pressure and temperature.
+
+        The behavior can be disabled by setting
+        ``params['initialize_operator_reference_from_initial_values'] = False``.
+
+        NOTE: This method is intended as a temporary bridge from PR #1696 until
+        downstream PRs on initialization have been merged.
+
+        """
+        if not self.params.get(
+            "initialize_operator_reference_from_initial_values", True
+        ):
+            return
+
+        for quantity_name in ("pressure", "temperature"):
+            variable_attr = f"{quantity_name}_variable"
+            if not hasattr(self, variable_attr):
+                continue
+
+            reference_value = cast(
+                pp.number, getattr(self.reference_variable_values, quantity_name, 0.0)
+            )
+            if np.isclose(reference_value, 0.0):
+                continue
+
+            variable_name = getattr(self, variable_attr)
+            domains = cast(list[pp.GridLike], self.mdg.subdomains())
+            variables = self.equation_system.get_variables([variable_name], domains)
+            if len(variables) == 0:
+                continue
+
+            values = self.equation_system.get_variable_values(
+                variables=variables, iterate_index=0
+            )
+            self.equation_system.set_variable_values(
+                values, variables=variables, reference=True
+            )
+
     def set_equation_system_manager(self) -> None:
         """Create an equation_system manager on the mixed-dimensional grid."""
         if not hasattr(self, "equation_system"):

src/porepy/numerics/ad/__init__.py

@@ -15,6 +15,7 @@
     equation_system,
     forward_mode,
     functions,
+    get_set_values,
     grid_operators,
     operator_functions,
     operators,
@@ -26,13 +27,15 @@
 from .equation_system import *
 from .forward_mode import *
 from .functions import *
+from .get_set_values import *
 from .grid_operators import *
 from .operator_functions import *
 from .operators import *
 from .surrogate_operator import *
 from .time_derivatives import *
 
 __all__.extend(ad_utils.__all__)
+__all__.extend(get_set_values.__all__)
 __all__.extend(operators.__all__)
 __all__.extend(operator_functions.__all__)
 __all__.extend(discretizations.__all__)

src/porepy/numerics/ad/_ad_parser.py

@@ -197,7 +197,7 @@ def _evaluate_single(
         #    them according to the operator.
         if op.is_leaf():
             if isinstance(op, pp.ad.MixedDimensionalVariable):
-                if op.is_previous_iterate or op.is_previous_time:
+                if op.is_previous_iterate or op.is_previous_time or op.is_reference:
                     # Empty vector like the global vector of unknowns for prev time/iter
                     # insert the values at the right dofs and slice.
                     vals = np.empty_like(
@@ -218,7 +218,7 @@ def _evaluate_single(
             # Atomic variables.
             elif isinstance(op, pp.ad.Variable):
                 # If a variable represents a previous iteration or time, parse values.
-                if op.is_previous_iterate or op.is_previous_time:
+                if op.is_previous_iterate or op.is_previous_time or op.is_reference:
                     return op.parse(equation_system.mdg)
                 # Otherwise use the current time and iteration values.
                 else: