Nonmatching grid tests #1334

zhangyh0713 · 2025-02-17T09:12:28Z

Proposed changes

Issue: #1312
We created a new class to generate non-matching grids, which can handle cases with 0 to 2 orthogonal fractures. It returns a refined fracture set and a refined interface map. Its functionality is tested in test_fluid_mass_balance.py. The fluid mass problem, previously tested only on matching grids, is also tested with non-matching grids.

Types of changes

What types of changes does this PR introduce to PorePy?
Put an x in the boxes that apply.

Minor change (e.g., dependency bumps, broken links).
Bugfix (non-breaking change which fixes an issue).
New feature (non-breaking change which adds functionality).
Breaking change (fix or feature that would cause existing functionality to not work as expected).
Testing (contribution related to testing of existing or new functionality).
Documentation (contribution related to adding, improving, or fixing documentation).
Maintenance (e.g., improve logic and performance, remove obsolete code).
Other:

Checklist

Put an x in the boxes that apply or explain briefly why the box is not relevant.

The documentation is up-to-date.
Static typing is included in the update.
This PR does not duplicate existing functionality.
The update is covered by the test suite (including tests added in the PR).
If new skipped tests have been introduced in this PR, pytest was run with the --run-skipped flag.

…d minor refactor

src/porepy/applications/md_grids/model_geometries.py

keileg

Nice work. I added some comments, and suggest to delete most of the documentation that reflect the development, rather than the result.

Main request: This function will be much more useful (also outside tests) if the user can controll the refinement process in fractures and interfaces. Please consider how this can be controlled by using parameters.

src/porepy/applications/md_grids/model_geometries.py

keileg · 2025-02-17T14:29:56Z

We should also make sure the geometry setup is directly tested, not just through the inclusion in higher-order tests.

IngridKJ · 2025-02-17T14:30:43Z

Note to selves after randomly meeting @keileg in the hallway: Make a test for testing the geometry itself.

…ractures class

keileg · 2025-02-20T13:06:08Z

This look fairly good now. Let's try to do a final change in today's session: We see if we can move the actual mesh generation to mdg_library, and let the class NonMatchingSquareDomainOrthhhogonalFractures fetch from the library.

Fixes to come

IvarStefansson · 2025-02-21T14:06:08Z

src/porepy/applications/md_grids/mdg_library.py

+        # For the fracture grids, we can use the pp.refinement.refine_grid_1d function
+        # and impose the ratio directly. For the interface grids, the simplest way of
+        # refining is to create a new mixed-dimensional grid with the same geometry, but
+        # with refined rates, and then map the new grids to the old ones.


What are "rates"?

IvarStefansson · 2025-02-21T14:12:19Z

src/porepy/applications/md_grids/model_geometries.py

+        # Refine and replace fracture grids. First we fetch the old fracture grids:
+        old_fracture_grids = self.mdg.subdomains(dim=1)
+
+        # Ratios which we want to refine the fracture grids with. Default, if no ratios
+        # are provided, is to refine by a factor of 2.
+        ratios = self.params.get(
+            "fracture_refinement_ratios", np.ones(len(old_fracture_grids) * 2)
+        )
+
+        # The actual refinement of the fracture grids.
+        new_fracture_grids = [
+            pp.refinement.refine_grid_1d(g=old_grid, ratio=ratio)
+            for old_grid, ratio in zip(old_fracture_grids, ratios)
+        ]
+
+        # Create a mapping between the old and new fracture grids.
+        grid_map = dict(zip(old_fracture_grids, new_fracture_grids))
+
+        # Refine and replace interface grids.
+        # Directly refining the already existing interface grids is not straightforward.
+        # Instead, we create a new mixed-dimensional grid and donate the interface grids
+        # to the original mixed-dimensional grid.
+
+        # First we gather the number of cells in x- and y- direction, as well as the
+        # physical dimensions of the original mixed-dimensional grid:
+        domain = pp.Domain(self.domain.bounding_box)
+        meshing_args = self.meshing_arguments()
+        nx_cells, phys_dims, _ = _preprocess_cartesian_args(domain, meshing_args, {})
+
+        # Then we create a new mdg with the same fractures and physical dimensions as
+        # the old one, but with a higher resolution:
+        nx_cells_new = np.array(nx_cells) * self.params.get(
+            "interface_refinement_ratio", 2
+        )
+        fracs = [self.fractures[i].pts for i in range(len(self.fractures))]
+        mdg_new = pp.meshing.cart_grid(fracs, nx_cells_new, physdims=phys_dims)
+
+        # Loop through all the interfaces in the new mixed-dimensional grid (donor) such
+        # that we can fill intf_map with a map between old and new interface grids.
+        intf_map: dict[
+            pp.MortarGrid,
+            Union[pp.MortarGrid, dict[mortar_grid.MortarSides, pp.Grid]],
+        ] = {}
+        for intf in mdg_new.interfaces(dim=1):
+            # Then, for each interface, we fetch the secondary grid which belongs to it.
+            _, g_sec = mdg_new.interface_to_subdomain_pair(intf)
+
+            # We then loop through all the interfaces in the original grid (recipient).
+            for intf_coarse in self.mdg.interfaces(dim=1):
+                # Fetch the secondary grid of the interface in the original grid.
+                _, g_sec_coarse = self.mdg.interface_to_subdomain_pair(intf_coarse)
+
+                # Checking the fracture number of the secondary grid in the recipient
+                # mdg. If they are the same, i.e., if the fractures are the same ones,
+                # we update the interface map.
+                if g_sec_coarse.frac_num == g_sec.frac_num:
+                    intf_map.update({intf_coarse: intf})
+
+        # Finally replace the subdomains and interfaces in the original
+        # mixed-dimensional grid.
+        self.mdg.replace_subdomains_and_interfaces(sd_map=grid_map, intf_map=intf_map)


Aren't lines 192 through 252 just a repetition of what happens during call to pp.mdg_library.square_with_orthogonal_fractures?

IvarStefansson · 2025-02-21T14:14:25Z

tests/applications/md_grids/test_mdg_library.py

+        # A matching grid would have equality here. We have refined the mortar grid,
+        # hence there should be more items in the projection matrix than there are
+        # fracture faces in the matrix grid.
+        assert non_zero_projection_primary > num_fracture_faces_from_matrix


For cartesian grids, we could check that these match the refinement ratios, right?

IvarStefansson · 2025-02-21T14:17:13Z

tests/applications/md_grids/test_model_geometries.py

+    params = {
+        "fracture_indices": fracture_indices,
+        "fracture_refinement_ratios": fracture_refinement_ratios,
+        "interface_refinement_ratio": interface_refinement_ratio,


Why is there a single interface ratio and multiple fracture ratios?

IvarStefansson

Some minor comments.

IngridKJ and others added 9 commits January 23, 2025 15:08

TST: Geometry mixin with nonmatching grids

6f50386

TST: Trying to make nonmatching interface grids (WIP)

156929e

TST: Example of creating nonmatching interface grid now works

bd5b8b3

DOC/ENH: Documentation in non-matching grid for applications/tests an…

6996974

…d minor refactor

MAINT: Refactor generation of mdg for application grid

a1dc1b6

TST: Preserve fracture number during grid refinement process

699a478

MAINT: Generatlity wrt number of fracs in nonmatching application grid

f75a0b1

TST: Run test_unit_conversion with nonmatching grid

e6066b4

MAINT: Remove code lines used while making nonmatching grid

2266ab3

zhangyh0713 requested review from keileg, IvarStefansson and jwboth as code owners February 17, 2025 09:12

keileg self-assigned this Feb 17, 2025