pmgbergen · Yuriyzabegaev · Nov 21, 2024 · Oct 24, 2024 · Nov 7, 2024 · Nov 8, 2024
@@ -5,11 +5,14 @@
 
 from __future__ import annotations
 
+import typing
+from pathlib import Path
 from typing import Optional
 
 import numpy as np
 
 import porepy as pp
+from porepy.fracs.fracture_network_2d import FractureNetwork2d
 
 
 def orthogonal_fractures_2d(
@@ -126,6 +129,35 @@ def benchmark_2d_case_3(size: pp.number = 1) -> list[pp.LineFracture]:
     return fractures
 
 
+def benchmark_2d_case_4() -> list[pp.LineFracture]:
+    """Returns a fracture set that corresponds to Case 4 from the 2D flow benchmark [1].
+
+    The setup is composed of 64 fractures grouped in 13 different connected networks,
+    ranging from isolated fractures up to tens of fractures each.
+
+    References:
+        - [1] Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N.,
+        Scotti, A., ... & Tatomir, A. (2018). Benchmarks for single-phase flow in
+        fractured porous media. Advances in Water Resources, 111, 239-258.
+
+    """
+    fracture_network_path = (
+        Path(__file__).parent
+        / "gmsh_file_library"
+        / "benchmark_2d_case_4"
+        / "fracture_network_benchmark_2d_case_4.csv"
+    )
+    # ``network_2d_from_csv`` is called with ``return_frac_id=False``, i.e.,
+    # ``fracture_network`` is actually a ``FractureNetwork2D`` object.
+    fracture_network = typing.cast(
+        FractureNetwork2d,
+        pp.fracture_importer.network_2d_from_csv(str(fracture_network_path)),
+    )
+    # We return only the fractures to stay consistent with the other functions' API from
+    # this file.
+    return fracture_network.fractures
+
+
 def seven_fractures_one_L_intersection(size: pp.number = 1) -> list[pp.LineFracture]:
     """Return a list of seven fractures with one L intersection.
 

@@ -0,0 +1,64 @@
+FID,START_X,START_Y,END_X,END_Y
+1,269.611206,152.05243,356.9240112,310.14123
+2,249.5117187,514.990780001,272.218872,470.97082
+3,258.3590698,515.574580001,271.9851684,490.9682
+4,270.6622924,524.702640001,269.1347046,147.78143
+5,355.8302002,348.479800001,337.5810733205,600
+6,366.9730835,338.132990001,426.9185141723,600
+7,198.237915,222.724420001,175.1561889,597.603030001
+8,151.2785034,261.724610001,154.4623059774,600
+9,29.5026855,300.724610001,96.3599853,514.82739
+10,386.0808105,33.3621800002,440.585083,275.191830001
+11,459.6350708,40.2413900001,461.751709,204.812620001
+12,297.180603,237.62103,468.1018066,40.2413900001
+13,312.5264892,272.01678,417.3016967,140.7832
+14,330.5181884,298.47522,439.5266723,156.6582
+15,340.5723877,320.70019,367.5598755,286.304380001
+16,492.9725952,312.762820001,576.5811157,419.6546
+17,505.6726684,309.05859,576.0520019,405.367190001
+18,537.4227905,297.94598,623.3187866,376.68463
+19,322.5338745,380.76941,521.8778076,593.552180001
+20,344.9320678,481.56122,409.8867798,503.959410001
+21,371.8098755,468.12219,510.6787109,383.009210001
+22,432.2849731,510.678830001,642.8280029,374.04999
+23,527.528634971,600,700,473.015615092
+24,0,333.73321,441.2443847,0
+25,13.4389038,342.692380001,347.171875,595.791990001
+26,22.3981933,450.203790001,311.3347778,291.176630001
+27,26.8778076,506.199220001,199.343811,400.92779
+28,44.7963867,528.597410001,365.0905151,342.692380001
+29,378.5294189,309.095210001,512.918518,116.470640001
+30,461.4027099,253.099610001,530.8370971,134.38922
+31,347.171875,374.04999,640.5881958,253.099610001
+32,490.5203857,268.77844,564.4343872,145.58844
+33,47.0361938,181.425410001,53.7556152,306.85541
+34,382.4152832,424.151000001,447.8997192,371.76343
+35,587.9967651,394.78222,549.1029663,362.635190001
+36,589.9812011,393.59161,527.6716919,313.8194
+37,597.125,378.90722,533.6248169,295.960200001
+38,533.6248169,448.75738,453.8527832,326.91638
+39,511.7966919,461.85419,489.5715942,395.17901
+40,565.3748779,425.34161,483.6184692,315.40698
+41,534.4185791,407.482240001,467.3466186,315.803830001
+42,627.2874756,527.3388,574.8999023,498.763610001
+43,644.3532104,519.00439,586.4093017,490.03241
+44,655.8626098,502.335630001,602.6812133,476.53863
+45,415.355896,585.679380001,391.9401855,561.47003
+46,417.3402099,578.535580001,397.8933105,554.326230001
+47,403.0526733,592.029420001,382.0183105,561.86682
+48,495.1278686,505.113580001,468.1403198,481.30121
+49,533.6248169,254.84381,420.9121093,159.196590001
+50,508.6217041,221.10943,441.152771,159.59363
+51,418.5308838,229.04681,312.961914,93.3154300004
+52,362.5714111,174.6748,322.883789,120.69983
+53,357.8088989,216.3468,295.102478,114.74658
+54,402.2589111,283.41882,366.1433105,226.66559
+55,337.5681762,253.256220001,374.4776001,211.18744
+56,386.7808838,264.765620001,509.8123169,101.25281
+57,473.2996826,278.65643,561.0092163,144.909240001
+58,471.7122192,253.653200001,554.6593017,129.034240001
+59,559.0249023,219.125,567.3593139,153.64044
+60,567.7561035,214.759400001,573.7092895,162.37182
+61,574.8999023,215.553040001,579.6624145,173.88104
+62,557.0404663,285.006410001,600.6968994,325.48761
+63,565.3748779,283.022030001,607.0468139,323.503230001
@@ -0,0 +1,86 @@
+"""
+This module contains the implementation of Case 4 from the 2D flow benchmark [1].
+
+The setup is composed of 64 fractures grouped in 13 different connected networks,
+ranging from isolated fractures up to tens of fractures each.
+
+Note:
+    At the current stage, the setup is meant only for performance profiling and does not
+    fully match the paper.
+
+References:
+    - [1] Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A.,
+      ... & Tatomir, A. (2018). Benchmarks for single-phase flow in fractured porous
+      media. Advances in Water Resources, 111, 239-258.
+
+"""
+
+import numpy as np
+
+import porepy as pp
+from porepy.models.protocol import PorePyModel
+
+solid_constants = pp.SolidConstants(residual_aperture=1e-4)
+
+
+class Geometry(PorePyModel):
+    """Geometry specification."""
+
+    def set_fractures(self) -> None:
+        """Setting a fracture list from the fracture set library."""
+        self._fractures = pp.applications.md_grids.fracture_sets.benchmark_2d_case_4()
+
+    @property
+    def domain(self) -> pp.Domain:
+        """Domain of the problem."""
+        # return pp.Domain({'xmax': 700, 'ymax': 600})
+        return pp.Domain({"xmax": 50, "ymax": 50})
+
+
+class BoundaryConditions(PorePyModel):
+    """Boundary conditions for Case 4 of the 2D flow benchmark.
+
+    Inflow on west (left) and prescribed pressure on east (right).
+
+    """
+
+    def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
+        """Pressure value of one atmosphere (101325 Pa) on west side."""
+        bounds = self.domain_boundary_sides(boundary_grid)
+        values = np.zeros(boundary_grid.num_cells)
+        values[bounds.west] = self.units.convert_units(101325, "Pa")
+        return values
+
+    def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
+        """Assign Dirichlet to the east and west boundary. The rest are Neumann by
+        default."""
+        bounds = self.domain_boundary_sides(sd)
+        bc = pp.BoundaryCondition(sd, bounds.east + bounds.west, "dir")
+        return bc
+
+    def bc_values_darcy_flux(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
+        """Inflow on the west boundary.
+
+        Per PorePy convention, the sign is negative for inflow and the value is
+        integrated over the boundary cell volumes. Since the inflow boundary contains a
+        fracture, the latter includes the fracture specific volume.
+
+        Parameters:
+            boundary_grid: Boundary grid.
+
+        Returns:
+            Boundary values.
+
+        """
+        bounds = self.domain_boundary_sides(boundary_grid)
+        values = np.zeros(boundary_grid.num_cells)
+        return values
+
+
+# Ignore type errors inherent to the ``SinglePhaseFlow`` class.
+class FlowBenchmark2dCase4Model(  # type: ignore[misc]
+    Geometry,
+    BoundaryConditions,
+    pp.fluid_mass_balance.SinglePhaseFlow,
+):
+    """Mixer class for case 4 from the 2d flow benchmark."""