Before rerun experiments

examples/.ipynb_checkpoints/mandel_model-checkpoint.py

@@ -0,0 +1,173 @@
+from dataclasses import dataclass
+from typing import Literal
+
+import numpy as np
+import porepy as pp
+import scipy.sparse
+from mandel_biot import MandelSetup as MandelSetup_
+from mandel_solvers import MandelSolverAssembler
+from porepy_common import PorepyNewtonSolver, PorepySimulation
+from scipy.sparse import csr_matrix
+
+from solver_selector.data_structures import ProblemContext
+from solver_selector.simulation_runner import Solver
+
+
+class MandelSetup(MandelSetup_):
+    def save_data_time_step(self) -> None:
+        # Prevent writing data that we don't need.
+        return pp.DataSavingMixin.save_data_time_step(self)
+
+
+@dataclass(frozen=True, kw_only=True, slots=True)
+class MandelContext(ProblemContext):
+    mat_size: int
+    time_step_value: float
+    cfl_max: float
+    cfl_mean: float
+
+    def get_array(self):
+        return np.array(
+            [
+                np.log(self.time_step_value),
+                np.log(self.mat_size),
+                np.log(self.cfl_max),
+                np.log(self.cfl_mean),
+            ],
+            dtype=float,
+        )
+
+
+class MandelSimulationModel(PorepySimulation):
+    def __init__(self, porepy_setup: MandelSetup):
+        self.porepy_setup: MandelSetup = porepy_setup
+        self.solver_assembler = MandelSolverAssembler(self)
+
+    def get_context(self) -> ProblemContext:
+        CFL_max, CFL_mean = self._compute_cfl()
+        return MandelContext(
+            mat_size=self.porepy_setup.equation_system._variable_num_dofs.sum(),
+            time_step_value=self.porepy_setup.time_manager.dt,
+            cfl_max=CFL_max,
+            cfl_mean=CFL_mean,
+        )
+
+    def assemble_solver(self, solver_config: dict) -> Solver:
+        linear_solver = self.solver_assembler.assemble_linear_solver(solver_config)
+        return PorepyNewtonSolver(linear_solver)
+
+    def get_primary_secondary_indices(
+        self, primary_variable: str
+    ) -> tuple[np.ndarray, np.ndarray]:
+        displacement = np.concatenate(
+            [
+                self.porepy_setup.equation_system.assembled_equation_indices[i]
+                for i in (
+                    "momentum_balance_equation",
+                    "normal_fracture_deformation_equation",
+                    "tangential_fracture_deformation_equation",
+                    "interface_force_balance_equation",
+                )
+            ]
+        )
+        pressure = np.concatenate(
+            [
+                self.porepy_setup.equation_system.assembled_equation_indices[i]
+                for i in (
+                    "mass_balance_equation",
+                    "interface_darcy_flux_equation",
+                    "well_flux_equation",
+                )
+            ]
+        )
+        if primary_variable == "displacement":
+            return displacement, pressure
+        elif primary_variable == "pressure":
+            return pressure, displacement
+        raise ValueError(f"{primary_variable=}")
+
+    def get_fixed_stress_stabilization(self, l_factor: float) -> csr_matrix:
+        porepy_setup = self.porepy_setup
+        mu_lame = porepy_setup.solid.shear_modulus()
+        lambda_lame = porepy_setup.solid.lame_lambda()
+        alpha_biot = porepy_setup.solid.biot_coefficient()
+        dim = 2
+
+        l_phys = alpha_biot**2 / (2 * mu_lame / dim + lambda_lame)
+        l_min = alpha_biot**2 / (4 * mu_lame + 2 * lambda_lame)
+
+        val = l_min * (l_phys / l_min) ** l_factor
+
+        diagonal_approx = val
+        subdomains = porepy_setup.mdg.subdomains()
+        cell_volumes = subdomains[0].cell_volumes
+        diagonal_approx *= cell_volumes
+
+        density = (
+            porepy_setup.fluid_density(subdomains)
+            .value(porepy_setup.equation_system)
+        )
+        diagonal_approx *= density
+
+        dt = porepy_setup.time_manager.dt
+        diagonal_approx /= dt
+
+        return scipy.sparse.diags(diagonal_approx)
+
+
+time_manager = pp.TimeManager(
+    schedule=[0, 1e3],
+    # schedule=[0, 1e2],
+    # schedule=[0, 10],
+    dt_init=10,
+    constant_dt=True,
+)
+
+units = pp.Units(
+    # m=1e-3,
+    # kg=1e3,
+)
+
+ls = 1 / units.m  # length scaling
+
+mandel_solid_constants = {
+    "lame_lambda": 1.65e9 * 1e-6,  # [MPa]
+    "shear_modulus": 2.475e9 * 1e-6,  # [MPa]
+    "specific_storage": 6.0606e-11 * 1e6,  # [MPa^-1]
+    "permeability": 9.869e-14 * 1e6,  # [mm^2]
+    "biot_coefficient": 1.0,  # [-]
+}
+
+mandel_fluid_constants = {
+    "density": 1e3,  # [kg * m^-3]
+    "viscosity": 1e-3,  # [Pa * s]
+}
+
+
+def make_mandel_setup(
+    model_size: Literal["small", "medium", "large"]
+) -> MandelSimulationModel:
+    if model_size == "small":
+        cell_size = 1
+    elif model_size == "medium":
+        cell_size = 0.75
+    elif model_size == "large":
+        cell_size = 0.5
+    else:
+        raise ValueError(model_size)
+
+    porepy_setup = MandelSetup(
+        {
+            "meshing_arguments": {
+                "cell_size": cell_size * ls,
+            },
+            "time_manager": time_manager,
+            "material_constants": {
+                "solid": pp.SolidConstants(mandel_solid_constants),
+                "fluid": pp.FluidConstants(mandel_fluid_constants),
+            },
+        }
+    )
+
+    porepy_setup.prepare_simulation()
+    return MandelSimulationModel(porepy_setup)

examples/.ipynb_checkpoints/mandel_solvers-checkpoint.py

@@ -0,0 +1,153 @@
+from typing import TYPE_CHECKING, Literal, Sequence
+
+from solvers_common import (
+    DirectSolverNode,
+    LinearSolver,
+    LinearSolverNames,
+    Preconditioner,
+    SolverAssembler,
+    SplittingSchur,
+)
+
+from solver_selector.solver_space import (
+    ConstantNode,
+    KrylovSolverNode,
+    NumericalParameter,
+    ParametersNode,
+    SolverConfigNode,
+    SplittingNode,
+)
+
+if TYPE_CHECKING:
+    from mandel_model import MandelSimulationModel
+
+
+#  ------------------- Solver space -------------------
+
+
+class SplittingNames:
+    fixed_stress = "splitting_fixed_stress"
+
+    primary_subsolver = "primary"
+    secondary_subsolver = "secondary"
+
+
+class FixedStressNode(SolverConfigNode):
+    type = SplittingNames.fixed_stress
+
+    def __init__(
+        self,
+        solver_nodes: dict[str, Sequence[SolverConfigNode]],
+        l_factor: float | NumericalParameter,
+    ) -> None:
+        self.l_factor = l_factor
+        self.solver_nodes = solver_nodes
+        splitting = SplittingNode.from_solvers_list(
+            solver_nodes=solver_nodes,
+            other_children=[
+                ParametersNode({"l_factor": self.l_factor}),
+                ConstantNode("primary_variable", "displacement"),
+                ConstantNode("method", "upper"),
+            ],
+        )
+        super().__init__(children=[splitting])
+
+    def copy(self):
+        return type(self)(
+            solver_nodes=self.solver_nodes,
+            l_factor=self.l_factor,
+        )
+
+
+def make_mandel_solver_space(l_factor: Literal["0", "1", "dynamic"]):
+    if l_factor == "0":
+        l_factor_val = 0
+    elif l_factor == "1":
+        l_factor_val = 1
+    elif l_factor == "dynamic":
+        l_factor_val = NumericalParameter(bounds=(0, 1), default=0, is_optimized=True)
+
+    prec_primary = DirectSolverNode()
+    prec_secondary = DirectSolverNode()
+
+    return KrylovSolverNode.from_preconditioners_list(
+        name=LinearSolverNames.gmres,
+        preconditioners=[
+            FixedStressNode(
+                {
+                    SplittingNames.primary_subsolver: [prec_primary],
+                    SplittingNames.secondary_subsolver: [prec_secondary],
+                },
+                l_factor=l_factor_val,
+            ),
+        ],
+        other_children=[ParametersNode({"tol": 1e-10, "restart": 30, "maxiter": 10})],
+    )
+
+
+#  ------------------- Solvers implementations -------------------
+
+
+class SplittingFixedStress(SplittingSchur):
+    def __init__(
+        self,
+        mandel_problem: "MandelSimulationModel",
+        linear_solver_primary: LinearSolver,
+        linear_solver_secondary: LinearSolver,
+        config: dict,
+    ):
+        self.simulation: "MandelSimulationModel" = mandel_problem
+        for param in ["primary_variable", "l_factor"]:
+            assert param in config
+        super().__init__(linear_solver_primary, linear_solver_secondary, config)
+
+    def update(self, mat):
+        super().update(mat)
+        prim_ind, sec_ind = self.simulation.get_primary_secondary_indices(
+            primary_variable=self.config["primary_variable"]
+        )
+        self.primary_ind = prim_ind
+        self.secondary_ind = sec_ind
+
+        mat_00 = mat[prim_ind[:, None], prim_ind]
+        mat_01 = mat[prim_ind[:, None], sec_ind]
+        mat_10 = mat[sec_ind[:, None], prim_ind]
+        mat_11 = mat[sec_ind[:, None], sec_ind]
+
+        l_factor = self.config["l_factor"]
+        stab = self.simulation.get_fixed_stress_stabilization(l_factor)
+        mat_S = mat_11 + stab
+
+        self.linear_solver_primary.update(mat_00)
+        self.linear_solver_secondary.update(mat_S)
+        self.mat_01 = mat_01
+        self.mat_10 = mat_10
+
+
+#  ------------------- Solver assembling routines -------------------
+
+
+class MandelSolverAssembler(SolverAssembler):
+    def assemble_preconditioner(
+        self,
+        preconditioner_config: dict,
+    ) -> Preconditioner:
+        preconditioner_type, preconditioner_params = list(
+            preconditioner_config.items()
+        )[0]
+        if preconditioner_type == SplittingNames.fixed_stress:
+            primary = SplittingNames.primary_subsolver
+            secondary = SplittingNames.secondary_subsolver
+            linear_solver_primary = self.assemble_linear_solver(
+                linear_solver_config=preconditioner_params[primary]
+            )
+            linear_solver_secondary = self.assemble_linear_solver(
+                linear_solver_config=preconditioner_params[secondary]
+            )
+            return SplittingFixedStress(
+                mandel_problem=self.simulation,
+                linear_solver_primary=linear_solver_primary,
+                linear_solver_secondary=linear_solver_secondary,
+                config=preconditioner_params,
+            )
+        return super().assemble_preconditioner(preconditioner_config)

examples/.ipynb_checkpoints/run-checkpoint.bash

@@ -0,0 +1,19 @@
+#!/bin/bash
+
+if [ $# -eq 0 ]; then
+    echo "Pass a number of how many times to execute each experiment. Example: \"bash run.bash 3\""
+    echo "Exiting."
+    exit 1
+fi
+
+cd 1
+# bash run.bash $1
+
+cd ../2_source_location
+bash run.bash $1
+
+cd ../3_exploration_sensitivity
+bash run.bash $1
+
+cd ../4_extended_solver_space
+# bash run.bash $1