Use proper mixins in concrete models

src/porepy/applications/boundary_conditions/model_boundary_conditions.py

@@ -9,7 +9,7 @@
 import porepy as pp
 
 
-class BoundaryConditionsMassDirWestEast(pp.BoundaryConditionMixin):
+class BoundaryConditionsMassDirWestEast(pp.PorePyModel):
     """Boundary conditions for the flow problem.
 
     Dirichlet boundary conditions are defined on the west and east boundaries. Some
@@ -73,7 +73,7 @@ def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.west + domain_sides.east, "dir")
 
 
-class BoundaryConditionsMassDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsMassDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the flow problem.
 
     Dirichlet boundary conditions are defined on the north and south boundaries. Some
@@ -137,7 +137,7 @@ def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.north + domain_sides.south, "dir")
 
 
-class BoundaryConditionsEnergyDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsEnergyDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the thermal problem.
 
     Dirichlet boundary conditions are defined on the north and south boundaries. Some
@@ -183,7 +183,7 @@ def bc_type_enthalpy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.north + domain_sides.south, "dir")
 
 
-class BoundaryConditionsMechanicsDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsMechanicsDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the mechanics with Dirichlet conditions on north and
     south boundaries.
 

src/porepy/applications/md_grids/model_geometries.py

@@ -7,7 +7,7 @@
 from . import domains, fracture_sets
 
 
-class SquareDomainOrthogonalFractures(pp.ModelGeometry):
+class SquareDomainOrthogonalFractures(pp.PorePyModel):
     """Create a mixed-dimensional grid for a square domain with up to two
     orthogonal fractures.
 
@@ -52,7 +52,7 @@ def set_domain(self) -> None:
         self._domain = domains.nd_cube_domain(2, self.domain_size)
 
 
-class CubeDomainOrthogonalFractures(pp.ModelGeometry):
+class CubeDomainOrthogonalFractures(pp.PorePyModel):
     """Create a mixed-dimensional grid for a cube domain with up to three
     orthogonal fractures.
 
@@ -78,7 +78,7 @@ def set_domain(self) -> None:
         self._domain = domains.nd_cube_domain(3, self.domain_size)
 
 
-class RectangularDomainThreeFractures(pp.ModelGeometry):
+class RectangularDomainThreeFractures(pp.PorePyModel):
     """A rectangular domain with up to three fractures.
 
     The domain is `[0, 2] x [0, 1]`.

src/porepy/applications/test_utils/models.py

@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 import inspect
-from typing import Any, Callable
+from typing import Any, Callable, cast
 
 import numpy as np
 
@@ -72,9 +72,7 @@ class Thermoporomechanics(  # type: ignore[misc]
     """Combine components needed for poromechanics simulation."""
 
 
-def model(
-    model_type: str, dim: int, num_fracs: int = 1
-) -> MassBalance | MomentumBalance | MassAndEnergyBalance | Poromechanics:
+def model(model_type: str, dim: int, num_fracs: int = 1) -> pp.PorePyModel:
     """Setup for tests."""
     # Suppress output for tests
     fracture_indices = [i for i in range(num_fracs)]
@@ -114,24 +112,18 @@ class Model(geometry, model_class):
         pass
 
     # Create an instance of the combined class
-    model = Model(params)
+    model = cast(pp.PorePyModel, Model(params))
 
     # Prepare the simulation
     # (create grids, variables, equations, discretize, etc.)
     model.prepare_simulation()
     return model
 
 
-class RobinDirichletNeumannConditions:
+class RobinDirichletNeumannConditions(pp.PorePyModel):
     """Mixin for applying Neumann, Dirichlet and Robin conditions for a
     thermoporomechanics model."""
 
-    params: dict
-
-    domain_boundary_sides: Callable[[pp.GridLike], pp.domain.DomainSides]
-
-    nd: int
-
     def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         """Assigns pressure values on the north and south boundary."""
         p_north = self.params.get("pressure_north", 1)
@@ -259,7 +251,7 @@ def subdomains_or_interfaces_from_method_name(
     return domains
 
 
-def _add_mixin(mixin, parent):
+def _add_mixin(mixin: type, parent: type) -> type:
     """Helper method to dynamically construct a class by adding a mixin.
 
     Multiple mixins can be added by nested calls to this method.
@@ -378,7 +370,7 @@ def compare_values(
         assert np.isclose(np.sum(values_0 - values_1), 0, atol=1e-10 + rtol)
 
 
-def get_model_methods_returning_ad_operator(model_setup) -> list[str]:
+def get_model_methods_returning_ad_operator(model_setup: pp.PorePyModel) -> list[str]:
     """Get all possible testable methods to be used in test_ad_operator_methods_xx.
 
     A testable method is one that:

src/porepy/applications/test_utils/well_models.py

@@ -1,18 +1,14 @@
 """Contains code for setting up a simple but non-trivial model with a well.
 """
 
+from typing import Literal
+
 import numpy as np
 
 import porepy as pp
 
 
-class OneVerticalWell:
-
-    domain: pp.Domain
-    """Domain for the model."""
-
-    units: pp.Units
-    """Simulation units provided by the solution strategy mixin."""
+class OneVerticalWell(pp.PorePyModel):
 
     def set_well_network(self) -> None:
         """Assign well network class."""
@@ -36,11 +32,11 @@ def meshing_arguments(self) -> dict:
 
         return mesh_sizes
 
-    def grid_type(self) -> str:
+    def grid_type(self) -> Literal["simplex"]:
         return "simplex"
 
 
-class BoundaryConditionsWellSetup(pp.BoundaryConditionMixin):
+class BoundaryConditionsWellSetup(pp.PorePyModel):
     """Boundary conditions for the well setup."""
 
     def _bc_type(self, sd: pp.Grid, well_cond: str) -> pp.BoundaryCondition:
@@ -158,6 +154,7 @@ def bc_type_fourier_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
 
 
 class WellPermeability(pp.constitutive_laws.CubicLawPermeability):
+
     def permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Permeability [m^2].
 

src/porepy/examples/flow_benchmark_2d_case_4.py

@@ -20,7 +20,6 @@
 import porepy as pp
 from porepy.examples.flow_benchmark_2d_case_1 import FractureSolidConstants
 from porepy.models.constitutive_laws import DimensionDependentPermeability
-from porepy.models.protocol import PorePyModel
 
 solid_constants = FractureSolidConstants(
     residual_aperture=1e-2,  # m
@@ -30,7 +29,7 @@
 )
 
 
-class Geometry(PorePyModel):
+class Geometry(pp.PorePyModel):
     """Geometry specification."""
 
     def set_fractures(self) -> None:
@@ -43,7 +42,7 @@ def domain(self) -> pp.Domain:
         return pp.Domain({"xmax": 700, "ymax": 600})
 
 
-class BoundaryConditions(PorePyModel):
+class BoundaryConditions(pp.PorePyModel):
     """Boundary conditions for Case 4 of the 2D flow benchmark.
 
     Inflow on west (left) and prescribed pressure on east (right).
@@ -114,6 +113,6 @@ class FlowBenchmark2dCase4Model(  # type: ignore[misc]
     Geometry,
     BoundaryConditions,
     Permeability,
-    pp.fluid_mass_balance.SinglePhaseFlow,
+    pp.SinglePhaseFlow,
 ):
     """Mixer class for case 4 from the 2d flow benchmark."""

src/porepy/examples/flow_benchmark_3d_case_3.py

@@ -36,12 +36,9 @@
 )
 
 
-class Geometry(pp.ModelGeometry):
+class Geometry(pp.PorePyModel):
     """Define Geometry as specified in Section 5.3 of the benchmark study [1]."""
 
-    params: dict
-    """User-defined model parameters."""
-
     def set_geometry(self) -> None:
         """Create mixed-dimensional grid and fracture network."""
 
@@ -71,6 +68,7 @@ def set_geometry(self) -> None:
 
 
 class IntersectionPermeability(Permeability):
+
     def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Constant intersection permeability.
 
@@ -90,7 +88,7 @@ def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator
         return self.isotropic_second_order_tensor(subdomains, permeability)
 
 
-class BoundaryConditions(pp.BoundaryConditionMixin):
+class BoundaryConditions(pp.PorePyModel):
     """Define inlet and outlet boundary conditions as specified by the benchmark."""
 
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:

src/porepy/examples/mandel_biot.py

@@ -33,13 +33,11 @@
 import scipy.optimize as opt
 
 import porepy as pp
-import porepy.models.fluid_mass_balance as mass
 import porepy.models.poromechanics as poromechanics
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.derived_models.biot import BiotPoromechanics
 from porepy.numerics.linalg.matrix_operations import sparse_array_to_row_col_data
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -117,7 +115,7 @@ class MandelSaveData:
     """Current simulation time."""
 
 
-class MandelDataSaving(VerificationDataSaving):
+class MandelDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -1230,12 +1228,9 @@ def _plot_consolidation_degree(self):
 
 
 # -----> Geometry
-class MandelGeometry(pp.ModelGeometry):
+class MandelGeometry(pp.PorePyModel):
     """Class for setting up the rectangular geometry."""
 
-    params: dict
-    """Simulation model parameters."""
-
     def set_domain(self) -> None:
         """Set the domain."""
         ls = self.units.convert_units(1, "m")  # length scaling
@@ -1255,15 +1250,11 @@ def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
 
 
 # -----> Boundary conditions
-class MandelBoundaryConditionsMechanicsTimeDependent(
-    pp.momentum_balance.BoundaryConditionsMomentumBalance,
-):
+class MandelBoundaryConditionsMechanicsTimeDependent(pp.PorePyModel):
+
     exact_sol: MandelExactSolution
     """Exact solution object."""
 
-    params: dict
-    """Parameter dictionary of the verification setup."""
-
     def vertical_load(self):
         """Retrieve and scale applied force.
 
@@ -1288,8 +1279,11 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
             Vectorial boundary condition representation.
 
         """
-        # Inherit bc from parent class. This sets all bc faces as Dirichlet.
-        bc = super().bc_type_mechanics(sd=sd)
+
+        # NOTE see BC for momentum balance
+        boundary_faces = self.domain_boundary_sides(sd).all_bf
+        bc = pp.BoundaryConditionVectorial(sd, boundary_faces, "dir")
+        bc.internal_to_dirichlet(sd)
 
         # Get boundary sides, retrieve data dict, and bc object
         sides = self.domain_boundary_sides(sd)
@@ -1323,7 +1317,8 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
             the North side of the domain.
 
         """
-        bc_vals = super().bc_values_displacement(boundary_grid)
+
+        bc_vals = np.zeros((self.nd, boundary_grid.num_cells)).ravel("F")
 
         sides = self.domain_boundary_sides(boundary_grid)
         # Cells of the boundary grid are faces of the parent subdomain.
@@ -1337,7 +1332,8 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         return bc_vals
 
 
-class MandelBoundaryConditionsSinglePhaseFlow(mass.BoundaryConditionsSinglePhaseFlow):
+class MandelBoundaryConditionsSinglePhaseFlow(pp.PorePyModel):
+
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Define boundary condition types for the Darcy flux.
 
@@ -1379,21 +1375,6 @@ class MandelSolutionStrategy(poromechanics.SolutionStrategyPoromechanics):
 
     """
 
-    def __init__(self, params: dict) -> None:
-        """Constructor of the class.
-
-        Parameters:
-            params: Parameters of the verification setup.
-
-        """
-        super().__init__(params)
-
-        self.exact_sol: MandelExactSolution
-        """Exact solution object."""
-
-        self.results: list[MandelSaveData] = []
-        """List of stored results from the verification."""
-
     def set_materials(self):
         """Set material parameters.