MAINT: Make fracture axes scale with domain sizes

Author: IvarStefansson

src/porepy/applications/md_grids/model_geometries.py

@@ -350,6 +350,10 @@ def fracture_params(self) -> dict:
             - major_axis_angles: Major axis angles of the fractures (default [0.0, 0.0])
             - num_points: Number of points to define each fracture (default [10, 10])
 
+        The fracture axes are scaled by the minimum of the domain sizes. For adjusting
+        the fracture centers, the user should override the property
+        :meth:`fracture_centers`.
+
         Returns:
             A dictionary with fracture parameters.
 
@@ -362,23 +366,27 @@ def fracture_params(self) -> dict:
             "dip_angles": np.array([np.pi / 2, np.pi / 2]),
             "major_axis_angles": np.array([0.0, 0.0]),
         }
-        # Update with user-provided parameters, converting units as needed.
-        units = ["-", "-", "m", "rad", "rad", "rad"]
-        user_params = self.params.get("fracture_params", {})
-        for key, unit in zip(default_params.keys(), units):
-            if key in user_params:
-                if unit == "m":
-                    default_params[key] = self.units.convert_units(
-                        user_params[key], "m"
-                    )
-                else:
-                    default_params[key] = user_params[key]
+
         if "fracture_minor_axes" not in default_params:
             default_params["fracture_minor_axes"] = default_params[
                 "fracture_major_axes"
             ]
         return default_params
 
+    @property
+    def fracture_minor_axes(self) -> np.ndarray:
+        params = self.fracture_params()
+        # Scale minor axes by the minimum domain size.
+        size = min(self.domain_sizes())
+        return params["fracture_minor_axes"] * size
+
+    @property
+    def fracture_major_axes(self) -> np.ndarray:
+        params = self.fracture_params()
+        # Scale major axes by the minimum domain size.
+        size = min(self.domain_sizes())
+        return params["fracture_major_axes"] * size
+
     @property
     def fracture_centers(self) -> tuple[np.ndarray, np.ndarray]:
         dx, dy, dz = self.domain_sizes()
@@ -387,16 +395,17 @@ def fracture_centers(self) -> tuple[np.ndarray, np.ndarray]:
         return center_1, center_2
 
     def set_fractures(self):
-        """Set the two elliptic fractures."""
+        """Set the elliptic fractures as defined in the fracture parameters and the
+        fracture centers method."""
         self._fractures = []
         params = self.fracture_params()
         for i in range(params["num_fractures"]):
             f = pp.create_elliptic_fracture(
                 center=self.fracture_centers[i],
                 strike_angle=params["strike_angles"][i],
                 dip_angle=params["dip_angles"][i],
-                major_axis=params["fracture_major_axes"][i],
-                minor_axis=params["fracture_minor_axes"][i],
+                major_axis=self.fracture_major_axes[i],
+                minor_axis=self.fracture_minor_axes[i],
                 major_axis_angle=params["major_axis_angles"][i],
                 num_points=params["num_points"][i],
             )

src/porepy/examples/geothermal_reservoir.py

@@ -272,9 +272,11 @@ def before_nonlinear_loop(self) -> None:
         """Hook called before the nonlinear solver loop.
 
         Rediscretize flow problem at the onset of second time interval, when injection
-        starts. The tailoring is done to recompute the Darcy flux discretization only at
-        this specific time, since the discretization remains constant during the rest of
-        the simulation. See also time dependency in :bc_type_darcy_flux.
+        starts and the BC type changes from Neumann to Dirichlet.
+
+        The tailoring is done to recompute the Darcy flux discretization only at this
+        specific time, since the discretization remains constant during the rest of the
+        simulation. See also time dependency in :bc_type_darcy_flux.
 
         NOTE: Strictly speaking, we could be even more conservative and ensure we don't
         rediscretize if the time step needs to be recomputed due to nonlinear solver
@@ -441,9 +443,11 @@ class ConstraintLineSearchNonlinearSolver(
         }
     )
     # Define domain sizes (x, y, z) and fracture size.
-    size = 1e3  # [m]
-    fracture_size = 1.5e2  # [m]
-    domain_sizes = np.array([1.0 * size, 1.0 * size, 1.0 * size])  # [m]
+    length_scale = 1e3  # [m]
+    fracture_size = 0.15  # [-], fraction of length_scale
+    domain_sizes = np.array(
+        [1.0 * length_scale, 1.0 * length_scale, 1.0 * length_scale]
+    )  # [m]
     # Define model parameters. See file example_params.py and the mixins used to define
     # :class:`GeothermalReservoirWellBCs` for other options.
     model_params = {
@@ -474,8 +478,8 @@ class ConstraintLineSearchNonlinearSolver(
         # Set geometry and meshing related parameters.
         "grid_type": "simplex",
         "meshing_arguments": {
-            "cell_size": size / 5.0,
-            "cell_size_fracture": fracture_size / 2.5,
+            "cell_size": length_scale / 5.0,
+            "cell_size_fracture": fracture_size * length_scale / 2.5,
         },
         "fracture_params": {  # Other options are available in the geometry mixin.
             "fracture_major_axes": np.array((fracture_size, fracture_size * 1.2)),