Extending fluid concept to fit CF & Standard formulations

src/porepy/__init__.py

@@ -176,18 +176,22 @@
 from porepy.models.boundary_condition import BoundaryConditionMixin
 from porepy.models.geometry import ModelGeometry
 from porepy.models.units import Units
-from porepy.models.material_constants import (
-    FluidConstants,
-    SolidConstants,
-    MaterialConstants,
-)
-
 
 from porepy.viz.data_saving_model_mixin import DataSavingMixin
 from porepy.viz.diagnostics_mixin import DiagnosticsMixin
 from porepy.models.solution_strategy import SolutionStrategy
 from porepy.models import constitutive_laws
 
+# composite subpackage
+from . import compositional
+from porepy.compositional.materials import (
+    FluidConstants,
+    SolidConstants,
+    MaterialConstants,
+)
+from porepy.compositional.base import Component, Phase, Fluid
+from porepy.compositional.compositional_mixins import CompositionalVariables, FluidMixin
+
 # "Primary" models
 from porepy.models import fluid_mass_balance, momentum_balance
 
@@ -233,6 +237,3 @@
 from porepy.applications.boundary_conditions import model_boundary_conditions
 from porepy.applications import test_utils
 from porepy import applications
-
-# composite subpackage
-from . import compositional

src/porepy/applications/boundary_conditions/model_boundary_conditions.py

@@ -3,7 +3,6 @@
 from __future__ import annotations
 
 import warnings
-from typing import Any
 
 import numpy as np
 
@@ -20,8 +19,6 @@ class BoundaryConditionsMassDirWestEast(pp.BoundaryConditionMixin):
 
     """
 
-    fluid: pp.FluidConstants
-
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Boundary condition type for Darcy flux.
 
@@ -54,7 +51,9 @@ def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         """
         domain_sides = self.domain_boundary_sides(boundary_grid)
         values = np.zeros(boundary_grid.num_cells)
-        values[domain_sides.west + domain_sides.east] = self.fluid.pressure()
+        values[domain_sides.west + domain_sides.east] = (
+            self.fluid.reference_component.pressure
+        )
         return values
 
     def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
@@ -84,8 +83,6 @@ class BoundaryConditionsMassDirNorthSouth(pp.BoundaryConditionMixin):
 
     """
 
-    fluid: pp.FluidConstants
-
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Boundary condition type for Darcy flux.
 
@@ -118,7 +115,9 @@ def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         """
         domain_sides = self.domain_boundary_sides(boundary_grid)
         values = np.zeros(boundary_grid.num_cells)
-        values[domain_sides.north + domain_sides.south] = self.fluid.pressure()
+        values[domain_sides.north + domain_sides.south] = (
+            self.fluid.reference_component.pressure
+        )
         return values
 
     def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
@@ -190,15 +189,6 @@ class BoundaryConditionsMechanicsDirNorthSouth(pp.BoundaryConditionMixin):
 
     """
 
-    params: dict[str, Any]
-    """Model parameters."""
-    solid: pp.SolidConstants
-    """Solid parameters."""
-    fluid: pp.FluidConstants
-    """Fluid parameters."""
-    nd: int
-    """Number of dimensions."""
-
     def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         """Boundary condition type for mechanics.
 
@@ -251,8 +241,8 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         u_s = np.tile(
             self.params.get("u_south", np.zeros(self.nd)), (boundary_grid.num_cells, 1)
         ).T
-        values[:, sides.north] = self.solid.convert_units(u_n, "m")[:, sides.north]
-        values[:, sides.south] = self.solid.convert_units(u_s, "m")[:, sides.south]
+        values[:, sides.north] = self.units.convert_units(u_n, "m")[:, sides.north]
+        values[:, sides.south] = self.units.convert_units(u_s, "m")[:, sides.south]
         return values.ravel("F")
 
 
@@ -266,8 +256,8 @@ class TimeDependentMechanicalBCsDirNorthSouth(BoundaryConditionsMechanicsDirNort
     def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         """Displacement values.
 
-        Initial value is u_y = self.solid.fracture_gap() +
-        self.solid.maximum_elastic_fracture_opening() at north boundary. Adding it on
+        Initial value is u_y = self.solid.fracture_gap +
+        self.solid.maximum_elastic_fracture_opening at north boundary. Adding it on
         the boundary ensures a stress-free initial state, as it compensates for those
         two values corresponding to zero traction contact according to the class
         :class:`~porepy.models.constitutive_laws.FractureGap`. For positive times,
@@ -287,8 +277,7 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         # Add fracture width on top if there is a fracture.
         if len(self.mdg.subdomains()) > 1:
             frac_val = (
-                self.solid.fracture_gap()
-                + self.solid.maximum_elastic_fracture_opening()
+                self.solid.fracture_gap + self.solid.maximum_elastic_fracture_opening
             )
         else:
             frac_val = 0

src/porepy/applications/md_grids/model_geometries.py

@@ -7,7 +7,7 @@
 from . import domains, fracture_sets
 
 
-class SquareDomainOrthogonalFractures:
+class SquareDomainOrthogonalFractures(pp.ModelGeometry):
     """Create a mixed-dimensional grid for a square domain with up to two
     orthogonal fractures.
 
@@ -16,21 +16,6 @@ class SquareDomainOrthogonalFractures:
 
     """
 
-    params: dict
-    """Parameters for the model geometry. Entries relevant for this mixin are:
-        - domain_size: The side length of the square domain.
-        - fracture_indices: List of indices of fractures to be included in the grid.
-
-    """
-    units: pp.Units
-    """Units for the model geometry."""
-    solid: pp.SolidConstants
-    """Solid constant object that takes care of scaling of solid-related quantities.
-    Normally, this is set by a mixin of instance
-    :class:`~porepy.models.solution_strategy.SolutionStrategy`.
-
-    """
-
     @property
     def domain_size(self) -> pp.number:
         """Return the side length of the square domain.
@@ -40,7 +25,7 @@ def domain_size(self) -> pp.number:
 
         """
         # Scale by length unit.
-        return self.solid.convert_units(self.params.get("domain_size", 1.0), "m")
+        return self.units.convert_units(self.params.get("domain_size", 1.0), "m")
 
     def set_fractures(self) -> None:
         """Assigns 0 to 2 fractures to the domain.
@@ -67,7 +52,7 @@ def set_domain(self) -> None:
         self._domain = domains.nd_cube_domain(2, self.domain_size)
 
 
-class CubeDomainOrthogonalFractures:
+class CubeDomainOrthogonalFractures(pp.ModelGeometry):
     """Create a mixed-dimensional grid for a cube domain with up to three
     orthogonal fractures.
 
@@ -76,26 +61,11 @@ class CubeDomainOrthogonalFractures:
 
     """
 
-    params: dict
-    """Parameters for the model geometry. Entries relevant for this mixin are:
-        - domain_size: The side length of the cube domain.
-        - fracture_indices: List of indices of fractures to be included in the grid.
-
-    """
-    units: pp.Units
-    """Units for the model geometry."""
-    solid: pp.SolidConstants
-    """Solid constant object that takes care of scaling of solid-related quantities.
-    Normally, this is set by a mixin of instance
-    :class:`~porepy.models.solution_strategy.SolutionStrategy`.
-
-    """
-
     @property
     def domain_size(self) -> pp.number:
         """Return the side length of the cube domain."""
         # Scale by length unit.
-        return self.solid.convert_units(self.params.get("domain_size", 1.0), "m")
+        return self.units.convert_units(self.params.get("domain_size", 1.0), "m")
 
     def set_fractures(self) -> None:
         """Assigns 0 to 3 fractures."""
@@ -120,12 +90,9 @@ class RectangularDomainThreeFractures(pp.ModelGeometry):
 
     """
 
-    params: dict
-    """Parameters for the model."""
-
     def set_fractures(self) -> None:
         # Length scale:
-        ls = self.solid.convert_units(1, "m")
+        ls = self.units.convert_units(1, "m")
 
         fracture_indices = self.params.get("fracture_indices", [0])
         fractures = [
@@ -137,7 +104,7 @@ def set_fractures(self) -> None:
 
     def meshing_arguments(self) -> dict:
         # Divide by length scale:
-        ls = self.solid.convert_units(1, "m")
+        ls = self.units.convert_units(1, "m")
 
         mesh_sizes = {
             # Cartesian: 2 by 8 cells.
@@ -157,7 +124,7 @@ def set_domain(self) -> None:
             self.params["grid_type"] = "cartesian"
 
         # Length scale:
-        ls = self.solid.convert_units(1, "m")
+        ls = self.units.convert_units(1, "m")
 
         # Mono-dimensional grid by default
         phys_dims = np.array([2, 1]) * ls
@@ -179,7 +146,7 @@ class OrthogonalFractures3d(CubeDomainOrthogonalFractures):
 
     def meshing_arguments(self) -> dict:
         # Length scale:
-        ls = self.solid.convert_units(1, "m")
+        ls = self.units.convert_units(1, "m")
 
         mesh_sizes = {
             "cell_size": 0.5 * ls,

src/porepy/applications/test_utils/models.py

@@ -15,7 +15,11 @@
 
 
 class NoPhysics(  # type: ignore[misc]
-    pp.ModelGeometry, pp.SolutionStrategy, pp.DataSavingMixin, pp.BoundaryConditionMixin
+    pp.ModelGeometry,
+    pp.SolutionStrategy,
+    pp.DataSavingMixin,
+    pp.BoundaryConditionMixin,
+    pp.FluidMixin,
 ):
     """A model with no physics, for testing purposes.
 
@@ -245,11 +249,12 @@ def subdomains_or_interfaces_from_method_name(
     assert len(signature.parameters) == 1
 
     # The domain is a list of either subdomains or interfaces.
+    domains: list[pp.Grid] | list[pp.MortarGrid]
     if "subdomains" in signature.parameters or "domains" in signature.parameters:
         # If relevant, filter out the domains that are not to be tested.
         domains = mdg.subdomains(dim=domain_dimension)
     elif "interfaces" in signature.parameters:
-        domains = mdg.interfaces(dim=domain_dimension)  # type: ignore[assignment]
+        domains = mdg.interfaces(dim=domain_dimension)
 
     return domains
 
@@ -328,8 +333,8 @@ def compare_scaled_primary_variables(
             variables=[var_name], time_step_index=0
         )
         # Convert back to SI units.
-        values_0 = setup_0.fluid.convert_units(scaled_values_0, var_unit, to_si=True)
-        values_1 = setup_1.fluid.convert_units(scaled_values_1, var_unit, to_si=True)
+        values_0 = setup_0.units.convert_units(scaled_values_0, var_unit, to_si=True)
+        values_1 = setup_1.units.convert_units(scaled_values_1, var_unit, to_si=True)
         compare_values(values_0, values_1, cell_wise=cell_wise)
 
 
@@ -369,7 +374,7 @@ def compare_scaled_model_quantities(
             )
             # Convert back to SI units.
             value = method(domains).value(setup.equation_system)
-            values.append(setup.fluid.convert_units(value, method_unit, to_si=True))
+            values.append(setup.units.convert_units(value, method_unit, to_si=True))
         compare_values(values[0], values[1], cell_wise=cell_wise)
 
 
@@ -447,4 +452,73 @@ def get_model_methods_returning_ad_operator(model_setup) -> list[str]:
         ):
             testable_methods.append(method)
 
-    return testable_methods
+    # Appending testable methods of the fluid
+    fluid_methods = [
+        method for method in dir(model_setup.fluid) if not method.startswith("_")
+    ]
+    testable_fluid_methods: list[str] = []
+    # The basic flow model has no energy-related
+    skip_methods = ["specific_enthalpy", "thermal_conductivity"]
+    for method in fluid_methods:
+        if method in skip_methods:
+            continue
+        # Get method in callable form
+        callable_method = getattr(model_setup.fluid, method)
+
+        # Retrieve method signature via inspect
+        try:
+            signature = inspect.signature(callable_method)
+        except TypeError:
+            continue
+
+        # Append method to the `testable_methods` list if the conditions are met
+        if (
+            len(signature.parameters) == 1
+            and (
+                "subdomains" in signature.parameters
+                or "interfaces" in signature.parameters
+                or "domains" in signature.parameters
+            )
+            and (
+                "pp.ad.Operator" in signature.return_annotation
+                or "pp.ad.DenseArray" in signature.return_annotation
+            )
+        ):
+            testable_fluid_methods.append(f"fluid.{method}")
+
+    # The fluid is represented by a reference phase, whose saturation is 1
+    # I.e., it's thermodynamic properties are equal to the fluid properties
+    phase_methods = [
+        method
+        for method in dir(model_setup.fluid.reference_phase)
+        if not method.startswith("_")
+    ]
+    testable_phase_methods: list[str] = []
+    for method in phase_methods:
+        if method in skip_methods:
+            continue
+        # Get method in callable form
+        callable_method = getattr(model_setup.fluid.reference_phase, method)
+
+        # Retrieve method signature via inspect
+        try:
+            signature = inspect.signature(callable_method)
+        except TypeError:
+            continue
+
+        # Append method to the `testable_methods` list if the conditions are met
+        if (
+            len(signature.parameters) == 1
+            and (
+                "subdomains" in signature.parameters
+                or "interfaces" in signature.parameters
+                or "domains" in signature.parameters
+            )
+            and (
+                "pp.ad.Operator" in signature.return_annotation
+                or "pp.ad.DenseArray" in signature.return_annotation
+            )
+        ):
+            testable_phase_methods.append(f"fluid.reference_phase.{method}")
+
+    return testable_methods + testable_fluid_methods + testable_phase_methods