Use proper mixins in concrete models

src/porepy/applications/test_utils/models.py

@@ -112,7 +112,7 @@ class Model(geometry, model_class):
         pass
 
     # Create an instance of the combined class
-    model: pp.PorePyModel = cast(pp.PorePyModel, Model(params))
+    model = cast(pp.PorePyModel, Model(params))
 
     # Prepare the simulation
     # (create grids, variables, equations, discretize, etc.)

src/porepy/examples/mandel_biot.py

@@ -38,7 +38,6 @@
 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
@@ -116,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]
@@ -1384,7 +1383,6 @@ def set_materials(self):
         """
         super().set_materials()
         self.exact_sol = MandelExactSolution(self)
-        self.results = []
 
         # Biot's coefficient must be one
         assert self.solid.biot_coefficient == 1

src/porepy/examples/terzaghi_biot.py

@@ -44,13 +44,10 @@
 import numpy as np
 
 import porepy as pp
-import porepy.models.fluid_mass_balance as mass
-import porepy.models.momentum_balance as mechanics
 import porepy.models.poromechanics as poromechanics
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.derived_models.biot import BiotPoromechanics
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -100,7 +97,7 @@ class TerzaghiSaveData:
     """Current simulation time."""
 
 
-class TerzaghiDataSaving(VerificationDataSaving):
+class TerzaghiDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     exact_sol: TerzaghiExactSolution
@@ -524,7 +521,7 @@ def meshing_arguments(self) -> dict:
 
 
 # -----> Boundary conditions
-class TerzaghiBoundaryConditionsMechanics(mechanics.BoundaryConditionsMomentumBalance):
+class TerzaghiBoundaryConditionsMechanics(pp.PorePyModel):
 
     def applied_load(self) -> pp.number:
         """Obtain vertical load in scaled [Pa]."""
@@ -542,8 +539,10 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
             the South, and rollers on the sides.
 
         """
-        # Inherit bc from parent class. This sets all bc faces as Dirichlet.
-        bc = super().bc_type_mechanics(sd=sd)
+        # Start with all faces as Dirichlet faces, analogous to base mechanics set-up.
+        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
         _, east, west, north, *_ = self.domain_boundary_sides(sd)
@@ -582,9 +581,8 @@ def bc_values_stress(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         return bc_values.ravel("F")
 
 
-class TerzaghiBoundaryConditionsFlow(
-    mass.BoundaryConditionsSinglePhaseFlow,
-):
+class TerzaghiBoundaryConditionsFlow(pp.PorePyModel):
+
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Define boundary condition types for the Darcy flux.
 
@@ -638,21 +636,6 @@ class TerzaghiSolutionStrategy(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: TerzaghiExactSolution
-        """Exact solution object."""
-
-        self.results: list[TerzaghiSaveData] = []
-        """List of stored results from the verification."""
-
     def set_materials(self):
         """Set material parameters.
 

src/porepy/models/constitutive_laws.py

@@ -3943,7 +3943,7 @@ def _mpsa_consistency(
         return consistency
 
 
-class BiotPoroMechanicsPorosity(PoroMechanicsPorosity):
+class BiotPoroMechanicsPorosity(pp.PorePyModel):
     """Porosity for poromechanical models following classical Biot's theory.
 
     The porosity is defined such that, after the chain rule is applied to the

src/porepy/models/derived_models/biot.py

@@ -77,35 +77,27 @@
 """
 
 import porepy as pp
-from porepy.models.poromechanics import (
-    ConstitutiveLawsPoromechanics,
+from porepy.models.poromechanics import Poromechanics
+
+
+class BiotPoromechanics(  # type: ignore[misc]
+    pp.constitutive_laws.SpecificStorage,
+    pp.constitutive_laws.BiotPoroMechanicsPorosity,
     Poromechanics,
-    SolutionStrategyPoromechanics,
-)
+):
+    """Biot-Poromechanics model with special constitutive laws for specific storage and
+    porosity.
 
+    The model is valid for single-phase, single-component and incompressible flow only.
+    A respective check is performed via overload of :meth:`set_materials`.
 
-class SolutionStrategyBiot(SolutionStrategyPoromechanics):
-    """Modified solution strategy for the Biot class"""
+    """
 
     def set_materials(self):
         """Set the material constants."""
         super().set_materials()
-        # Check that fluid compressibility is zero, otherwise Biot class doesn't hold
+        # Check that fluid compressibility is zero, otherwise Biot class doesn't hold.
         # NOTE Biot is tested for 1 phase 1 component.
         assert self.fluid.num_components == 1
         assert self.fluid.num_phases == 1
         assert self.fluid.reference_component.compressibility == 0
-
-
-class ConstitutiveLawsBiot(
-    pp.constitutive_laws.SpecificStorage,
-    pp.constitutive_laws.BiotPoroMechanicsPorosity,
-    ConstitutiveLawsPoromechanics,
-): ...
-
-
-class BiotPoromechanics(  # type: ignore[misc]
-    ConstitutiveLawsBiot,
-    SolutionStrategyBiot,
-    Poromechanics,
-): ...

src/porepy/models/protocol.py

@@ -16,7 +16,7 @@
 """
 
 from pathlib import Path
-from typing import TYPE_CHECKING, Callable, Literal, Optional, Protocol, Sequence
+from typing import TYPE_CHECKING, Any, Callable, Literal, Optional, Protocol, Sequence
 
 import numpy as np
 import scipy.sparse as sps
@@ -537,6 +537,9 @@ class SolutionStrategyProtocol(Protocol):
         """Time step as an automatic differentiation scalar."""
         nonlinear_solver_statistics: pp.SolverStatistics
         """Solver statistics for the nonlinear solver."""
+        results: list[Any]
+        """A list of results collected by the data saving mixin in
+        :meth:`~porepy.viz.data_saving_model_mixin.DataSavingMixin.collect_data`."""
 
         @property
         def time_step_indices(self) -> np.ndarray:

src/porepy/models/solution_strategy.py

@@ -7,7 +7,6 @@
 
 from __future__ import annotations
 
-import abc
 import logging
 import time
 import warnings
@@ -22,7 +21,7 @@
 logger = logging.getLogger(__name__)
 
 
-class SolutionStrategy(abc.ABC, pp.PorePyModel):
+class SolutionStrategy(pp.PorePyModel):
     """This is a class that specifies methods that a model must implement to
     be compatible with the linearization and time stepping methods.
 
@@ -114,6 +113,9 @@ def __init__(self, params: Optional[dict] = None):
         """Restart options. The template is provided in `SolutionStrategy.__init__`."""
         self.ad_time_step = pp.ad.Scalar(self.time_manager.dt)
         """Time step as an automatic differentiation scalar."""
+        self.results: list[Any] = []
+        """A list of results collected by the data saving mixin in
+        :meth:`~porepy.viz.data_saving_model_mixin.DataSavingMixin.collect_data`."""
 
         self.set_solver_statistics()
 

src/porepy/viz/data_saving_model_mixin.py

@@ -10,7 +10,7 @@
 from __future__ import annotations
 
 from pathlib import Path
-from typing import Optional, Union
+from typing import Any, Optional, Union
 
 import numpy as np
 
@@ -29,14 +29,19 @@ class DataSavingMixin(pp.PorePyModel):
     def save_data_time_step(self) -> None:
         """Export the model state at a given time step and log time.
 
-        The options for exporting times are:
-            * `None`: All time steps are exported
-            * `list`: Export if time is in the list. If the list is empty, then no
-            times are exported.
+        The options for exporting times can be given as ``params['times_to_export']``:
+
+        - ``None``: All time steps are exported.
+        - ``list``: Export if time is in the list. If the list is empty, then no
+          times are exported.
 
         In addition, save the solver statistics to file if the option is set.
 
+        Finally, :meth:`collect_data` is called and stored in :attr:`results` for
+        data collection and verification in runtime.
+
         """
+
         # Fetching the desired times to export.
         times_to_export = self.params.get("times_to_export", None)
         if times_to_export is None:
@@ -55,6 +60,40 @@ def save_data_time_step(self) -> None:
         # Save solver statistics to file.
         self.nonlinear_solver_statistics.save()
 
+        # Collecting and storing data in runtime for analysis. If default value of None
+        # is returned, nothing is stored to not burden memory.
+        if not self._is_time_dependent():  # stationary problem
+            if (
+                self.nonlinear_solver_statistics.num_iteration > 0
+            ):  # avoid saving initial condition
+                collected_data = self.collect_data()
+                if collected_data is not None:
+                    self.results.append(collected_data)
+        else:  # time-dependent problem
+            t = self.time_manager.time  # current time
+            scheduled = self.time_manager.schedule[1:]  # scheduled times except t_init
+            if any(np.isclose(t, scheduled)):
+                collected_data = self.collect_data()
+                if collected_data is not None:
+                    self.results.append(collected_data)
+
+    def collect_data(self) -> Any:
+        """Collect relevant simulation data to be stored in attr:`results`.
+
+        Override to collect data respectively. By default, this method returns None and
+        nothing is stored.
+
+        For stationary problems, this method is called in every iteration. For time
+        dependent problems, it is called after convergence of a time step which is
+        scheduled by the time manager.
+
+        Returns:
+            Any data structure relevant for future verification. By default, None.
+            If it is not None, it is stored in :attr:`results`.
+
+        """
+        return None
+
     def write_pvd_and_vtu(self) -> None:
         """Helper function for writing the .vtu and .pvd files and time information."""
         self.exporter.write_vtu(self.data_to_export(), time_dependent=True)
@@ -237,29 +276,3 @@ def load_data_from_pvd(
         self.time_manager.load_time_information(times_file)
         self.time_manager.set_time_and_dt_from_exported_steps(time_index)
         self.exporter._time_step_counter = time_index
-
-
-class VerificationDataSaving(DataSavingMixin):
-    """Class to store relevant data for a generic verification setup."""
-
-    results: list
-    """List of objects containing the results of the verification."""
-
-    def save_data_time_step(self) -> None:
-        """Save data to the `results` list."""
-        if not self._is_time_dependent():  # stationary problem
-            if (
-                self.nonlinear_solver_statistics.num_iteration > 0
-            ):  # avoid saving initial condition
-                collected_data = self.collect_data()
-                self.results.append(collected_data)
-        else:  # time-dependent problem
-            t = self.time_manager.time  # current time
-            scheduled = self.time_manager.schedule[1:]  # scheduled times except t_init
-            if any(np.isclose(t, scheduled)):
-                collected_data = self.collect_data()
-                self.results.append(collected_data)
-
-    def collect_data(self):
-        """Collect relevant data for the verification setup."""
-        raise NotImplementedError()

tests/applications/test_convergence_analysis.py

@@ -26,7 +26,6 @@
 )
 from porepy.models.fluid_mass_balance import SinglePhaseFlow
 from porepy.utils.txt_io import read_data_from_txt
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 
 # -----> Fixtures that are required on a module level.
@@ -461,7 +460,7 @@ class StationaryModelSaveData:
         error_var_0: float  # error associated with variable 0
         error_var_1: float  # error associated with variable 1
 
-    class StationaryModelDataSaving(VerificationDataSaving):
+    class StationaryModelDataSaving(pp.PorePyModel):
         """Class that collects and store data."""
 
         def collect_data(self) -> StationaryModelSaveData:
@@ -485,9 +484,7 @@ def collect_data(self) -> StationaryModelSaveData:
     class StationaryModelSolutionStrategy(pp.SolutionStrategy):
         """Solution strategy for the stationary flow model."""
 
-        def __init__(self, params: dict):
-            super().__init__(params)
-            self.results: list[StationaryModelSaveData] = []
+        results: list[StationaryModelSaveData]
 
         def _is_nonlinear_problem(self) -> bool:
             """Whether the model is non-linear."""
@@ -528,7 +525,7 @@ class TimeDependentModelSaveData:
         error_var_0: float  # error associated with variable 0
         error_var_1: float  # error associated with variable 1
 
-    class TimeDependentModelDataSaving(VerificationDataSaving):
+    class TimeDependentModelDataSaving(pp.PorePyModel):
         """Class that collects and store data."""
 
         def collect_data(self) -> TimeDependentModelSaveData:
@@ -552,9 +549,7 @@ def collect_data(self) -> TimeDependentModelSaveData:
     class TimeDependentModelSolutionStrategy(pp.SolutionStrategy):
         """Solution strategy for the time-dependent flow model."""
 
-        def __init__(self, params: dict):
-            super().__init__(params)
-            self.results: list[TimeDependentModelSaveData] = []
+        results: list[TimeDependentModelSaveData]
 
         def _is_nonlinear_problem(self) -> bool:
             """Whether the problem is non-linear."""

tests/functional/setups/manu_flow_comp_2d_frac.py

@@ -34,7 +34,6 @@
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 from tests.functional.setups.manu_flow_incomp_frac_2d import (
     ManuIncompSaveData,
     ManuIncompUtils,
@@ -75,7 +74,7 @@ class ManuCompSaveData(ManuIncompSaveData):
     """Current simulation time."""
 
 
-class ManuCompDataSaving(VerificationDataSaving):
+class ManuCompDataSaving(pp.PorePyModel):
     """Mixin class to store relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -703,12 +702,6 @@ def __init__(self, params: dict):
         """Constructor of the class."""
         super().__init__(params)
 
-        self.exact_sol: ManuCompExactSolution2d
-        """Exact solution object."""
-
-        self.results: list[ManuCompSaveData] = []
-        """Object that stores exact and approximated solutions and L2 errors."""
-
         self.subdomain_darcy_flux_variable: str = "darcy_flux"
         """Keyword to access the subdomain Darcy fluxes."""