ENH: nested grid clean up (#431)

Co-authored-by: Lars Hustoft (SUB CCI UKI) <lahus@st-lintgx0472.st.statoil.no>

Author: lahus-0

docs/nestedhybridgrid.rst

@@ -41,6 +41,10 @@ The example here runs within RMS, but similar workflows can be created for file
     # Load grid and region property
     grid = xtgeo.grid_from_roxar(project, "Simgrid")
     region = xtgeo.gridproperty_from_roxar(project, "Simgrid", "REGION")
+    
+    # Optionally load any other parameter you wish to preserve on grid
+    #zone = xtgeo.gridproperty_from_roxar(project, "Simgrid", "Zone")
+    #grid.append_prop(zone)
 
     # Create nested hybrid grid (e.g. refine region 2 by 2×2×1)
     merged, nnc_table = create_nested_hybrid_grid(
@@ -50,11 +54,12 @@ The example here runs within RMS, but similar workflows can be created for file
     # store merged grid in RMS (or file)
     merged.to_roxar(project, "NestedHybrid")
 
-    merged.get_prop_by_name("NEST_ID")
-    nest_id.to_roxar(project, "NestedHybrid", nest_id.name)
+    # Optionally extract and store any necessary parameters from grid
+    region2 = merged.get_prop_by_name("REGION")
+    region2.to_roxar(project, "NestedHybrid", region2.name)
 
     # write the NNC pandas to disk; this will be applied for computing NNC's in the next script
-    nnc_table.write_csv("path_to_some_csv_file.csv", index=False)
+    nnc_table.to_csv("path_to_some_csv_file.csv", index=False)
 
 
 The next step is to do a rescaling from the original geogrid to the merged grid
@@ -64,6 +69,7 @@ Further, we need to create NNC transmissibilities and generate file for flow sim
 
 .. code-block:: python
 
+    import pandas as pd
     import xtgeo
     from fmu.tools.nestedhybridgrid import (
         nnc_to_flowsimulator_input,
@@ -74,7 +80,6 @@ Further, we need to create NNC transmissibilities and generate file for flow sim
 
     # Load grid and region property which may be stored in RMS
     nested = xtgeo.grid_from_roxar(project, GNAME)
-    region = nested.get_prop_by_name("NEST_ID")  # if needed for QC
 
     # load the NNC table
     nnc_table = pd.read_csv("path_to_some_nnc_file.csv")
@@ -88,15 +93,15 @@ Further, we need to create NNC transmissibilities and generate file for flow sim
 
     # compute transmissibilities. Note that flow simulators do this for the normal cells/faults
     # so strictly speaking, only nnc_hybrid is needed here.
-    tranx, trany, tranz, nnc_fault, nnc_hybrid, rbnd = merged.get_transmissibilities(
+    tranx, trany, tranz, nnc_fault, nnc_hybrid, rbnd = nested.get_transmissibilities(
         permx, permy, permz, ntg, nnc_table=nnc_table
     )
 
     # Export NNC keyword for Eclipse / OPM Flow
     nnc_to_flowsimulator_input(nnc_hybrid, "some_path/NNC_HYBRID.INC")
 
     # Or map NNCs onto grid properties for visualisation
-    tx_nnc, ty_nnc, tz_nnc = nnc_to_gridproperty(merged, nnc_hybrid)
+    tx_nnc, ty_nnc, tz_nnc = nnc_to_gridproperty(nested, nnc_hybrid)
     tx_nnc.to_roxar(project, GNAME, "TRANX_NNC_QC")  # etc
 
 Concepts
@@ -127,15 +132,6 @@ This table is passed to :meth:`xtgeo.Grid.get_transmissibilities` via the
 face-overlap calculations (Sutherland–Hodgman algorithm) and two-point flux
 approximation (TPFA).
 
-NEST_ID property
-^^^^^^^^^^^^^^^^
-
-The merged grid carries a discrete property called ``NEST_ID``:
-
-- **0** — inactive hole cells (carved out to make room for the refined region)
-- **1** — mother (coarse) cells
-- **2** — refined cells
-
 Eclipse / OPM Flow export
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 

pyproject.toml

@@ -44,7 +44,7 @@ dependencies = [
     "pyyaml>=5.3",
     "scipy>=1.2",
     "xlrd>=1.2",
-    "xtgeo>=4.21.0",
+    "xtgeo>=4.23.0",
 ]
 
 [project.optional-dependencies]

src/fmu/tools/nestedhybridgrid/nestedhybrid.py

@@ -130,6 +130,8 @@ def _compute_nnc_table(
     crop_origin: tuple[int, int, int],
     refinement: tuple[int, int, int],
     coarse_ncol: int,
+    lmap1: np.ndarray,
+    lmap2: np.ndarray,
 ) -> pd.DataFrame:
     """Compute NNC cell-pair mapping between mother and refined cells.
 
@@ -155,6 +157,9 @@ def _compute_nnc_table(
         crop_origin: 0-based ``(i0, j0, k0)`` origin of the crop box.
         refinement: ``(rcol, rrow, rlay)`` refinement factors.
         coarse_ncol: Number of columns in the coarse grid (grid1 in the merge).
+        lmap1: Numpy array with layer_mapping (input k -> output k) for grid1
+        lmap2: Numpy array with layer_mapping (input k -> output k) for grid2
+
 
     Returns:
         A DataFrame with columns ``I1, J1, K1, I2, J2, K2, DIRECTION``.
@@ -163,6 +168,7 @@ def _compute_nnc_table(
 
     rcol, rrow, rlay = refinement
     i0, j0, k0 = crop_origin
+
     # In the merged grid, grid2 (refined) starts after a 1-column gap:
     i_offset = coarse_ncol + 1
 
@@ -234,10 +240,10 @@ def _compute_nnc_table(
                         {
                             "I1": mi + 1,
                             "J1": mj + 1,
-                            "K1": mk + 1,
+                            "K1": lmap1[mk] + 1,
                             "I2": ri + i_offset + 1,
                             "J2": rj + 1,
-                            "K2": rk + 1,
+                            "K2": lmap2[rk] + 1,
                             "DIRECTION": direction,
                         }
                     )
@@ -271,6 +277,38 @@ def _set_actnum_by_region(
     grid.set_actnum(actnum)
 
 
+def _generate_layer_mappings(
+    coarse_nlay: int,
+    refined_nlay: int,
+    refinement: tuple[int, int, int],
+    crop_origin: tuple[int, int, int],
+) -> tuple[np.ndarray, np.ndarray]:
+    """Generate mappings from old to new layer number.
+    Args:
+        coarse_nlay: Number of layers in the coarse grid (grid1 in the merge).
+        refined_nlay: Number of layers in the refined grid (grid2 in the merge).
+        crop_origin: 0-based ``(i0, j0, k0)`` origin of the crop box.
+        refinement: ``(rcol, rrow, rlay)`` refinement factors.
+
+    Returns:
+        lmap1: Numpy array with layer_mapping (input k -> output k) for grid1
+        lmap2: Numpy array with layer_mapping (input k -> output k) for grid2
+    """
+
+    _, _, rlay = refinement
+    _, _, k0 = crop_origin
+
+    lmap1 = np.arange(coarse_nlay, dtype=np.int32)
+    lmap1 = lmap1 + np.where(
+        lmap1 < k0,
+        0,
+        (rlay - 1) * np.minimum(int(refined_nlay / rlay), lmap1 - k0),
+    )
+    lmap2 = np.arange(refined_nlay, dtype=np.int32) + k0
+
+    return (lmap1, lmap2)
+
+
 # ---------------------------------------------------------------------------
 # Public API
 # ---------------------------------------------------------------------------
@@ -281,17 +319,16 @@ def create_nested_hybrid_grid(
     region: xtgeo.GridProperty,
     target_region_id: int,
     refinement: tuple[int, int, int],
-) -> tuple[xtgeo.Grid, pd.DataFrame]:
+) -> tuple[
+    xtgeo.Grid,
+    pd.DataFrame,
+]:
     """Create a nested hybrid grid by refining one region and merging it back.
 
     The cells belonging to *target_region_id* are replaced by a refined
-    (subdivided) version of the same region.  A ``NEST_ID`` discrete property
-    is attached to the merged grid, encoding the nested hybrid structure:
-
-    - ``NEST_ID == 1``: coarse (mother) grid cells.
-    - ``NEST_ID == 2``: refined grid cells.
+    (subdivided) version of the same region.
 
-    In addition, a **NNC mapping table** is returned that lists every
+    A **NNC mapping table** is returned that lists every
     mother ↔ refined cell pair that should be connected by a Non-Neighbour
     Connection (NNC).  The table is derived from the topological knowledge
     available at merge time (which original cell was refined and how its
@@ -316,9 +353,9 @@ def create_nested_hybrid_grid(
         refinement: ``(ncol, nrow, nlay)`` refinement factors.
 
     Returns:
-        A tuple ``(merged_grid, nnc_table)`` where *merged_grid* is a new
-        :class:`xtgeo.Grid` with the refined region stitched back into the
-        coarse grid and *nnc_table* is a :class:`pandas.DataFrame` mapping
+        A tuple ``(merged_grid, nnc_table)`` where *merged_grid*
+        is a new :class:`xtgeo.Grid` with the refined region stitched back into
+        the coarse grid and *nnc_table* is a :class:`pandas.DataFrame` mapping
         mother cells to their connected refined cells.
     """
     if any(r < 1 for r in refinement):
@@ -343,10 +380,19 @@ def create_nested_hybrid_grid(
     # 2. Refine the cropped grid.
     refined = cropped.copy()
     rcol, rrow, rlay = refinement
+    _, _, olay = crop_origin
     refined.refine(refine_col=rcol, refine_row=rrow, refine_layer=rlay)
     _logger.info("Refined cropped grid dimensions: %s", refined.dimensions)
 
-    # 3. Compute the NNC mapping table *before* deactivation mutates anything.
+    # 3. Generate layer mappings
+    lmap1, lmap2 = _generate_layer_mappings(
+        coarse_nlay=grid.nlay,
+        refined_nlay=refined.nlay,
+        crop_origin=crop_origin,
+        refinement=refinement,
+    )
+
+    # 4. Compute the NNC mapping table *before* deactivation mutates anything.
     #    This uses the original region property to find boundary faces and
     #    maps them through the crop → refine → merge index chain.
     nnc_table = _compute_nnc_table(
@@ -355,37 +401,20 @@ def create_nested_hybrid_grid(
         crop_origin=crop_origin,
         refinement=refinement,
         coarse_ncol=grid.ncol,
+        lmap1=lmap1,
+        lmap2=lmap2,
     )
 
-    # 4. Deactivate the target region in the coarse grid (will be replaced).
+    # 5. Deactivate the target region in the coarse grid (will be replaced).
     coarse_region = grid.get_prop_by_name(region.name)
     _set_actnum_by_region(grid, coarse_region, target_region_id, invert=False)
 
-    # 5. In the refined grid keep only target-region cells active.
+    # 6. In the refined grid keep only target-region cells active.
     refined_region = refined.get_prop_by_name(region.name)
     _set_actnum_by_region(refined, refined_region, target_region_id, invert=True)
 
-    # 6. Create NEST_ID properties before merging (1=mother, 2=refined).
-    nest_id_coarse = xtgeo.GridProperty(
-        grid,
-        name="NEST_ID",
-        discrete=True,
-        values=np.where(grid.get_actnum().values == 1, 1, 0).astype(np.int32),
-        codes={0: "inactive", 1: "mother", 2: "refined"},
-    )
-    grid.append_prop(nest_id_coarse)
-
-    nest_id_refined = xtgeo.GridProperty(
-        refined,
-        name="NEST_ID",
-        discrete=True,
-        values=np.where(refined.get_actnum().values == 1, 2, 0).astype(np.int32),
-        codes={0: "inactive", 1: "mother", 2: "refined"},
-    )
-    refined.append_prop(nest_id_refined)
-
     # 7. Merge the two grids.
-    merged = xtgeo.grid_merge(grid, refined)
+    merged = xtgeo.grid_merge(grid, refined, lmap1, lmap2)
     _logger.info("Merged grid dimensions: %s", merged.dimensions)
 
     return merged, nnc_table

tests/nestedhybridgrid/test_nestedhybrid.py

@@ -10,6 +10,9 @@
     nnc_to_flowsimulator_input,
     nnc_to_gridproperty,
 )
+from fmu.tools.nestedhybridgrid.nestedhybrid import (
+    _generate_layer_mappings,
+)
 
 # ---------------------------------------------------------------------------
 # Helpers
@@ -77,23 +80,23 @@ def test_basic_merge_dimensions(self):
         assert merged.nlay >= grid.nlay
 
     def test_nest_id_property_attached(self):
-        """The merged grid must have a NEST_ID property."""
+        """The merged grid must have a refinement region property."""
         grid, region, rid = _make_box_grid_with_region(dimension=(6, 6, 2))
         merged, _ = create_nested_hybrid_grid(grid, region, rid, refinement=(1, 1, 1))
 
-        nest_id = merged.get_prop_by_name("NEST_ID")
+        nest_id = merged.get_prop_by_name(region.name)
         assert nest_id is not None
         unique_vals = set(np.unique(np.ma.filled(nest_id.values, fill_value=0)))
         # Must contain at least mother (1) and refined (2) cells
         assert 1 in unique_vals
         assert 2 in unique_vals
 
     def test_nest_id_values_consistent(self):
-        """Active cells should only have NEST_ID in {1, 2}."""
+        """Active cells should only have refinement region in {1, 2}."""
         grid, region, rid = _make_box_grid_with_region(dimension=(6, 6, 2))
         merged, _ = create_nested_hybrid_grid(grid, region, rid, refinement=(1, 1, 1))
 
-        nest_id = merged.get_prop_by_name("NEST_ID")
+        nest_id = merged.get_prop_by_name(region.name)
         actnum = merged.get_actnum()
 
         active_mask = actnum.values == 1
@@ -125,20 +128,84 @@ def test_original_grid_not_mutated(self):
         orig_nactive = grid.nactive
         orig_region_sum = int(np.ma.filled(region.values, 0).sum())
 
-        _ = create_nested_hybrid_grid(grid, region, rid, refinement=(2, 2, 1))
+        _, _ = create_nested_hybrid_grid(grid, region, rid, refinement=(2, 2, 1))
 
         assert grid.ncol == orig_ncol
         assert grid.nactive == orig_nactive
         assert int(np.ma.filled(region.values, 0).sum()) == orig_region_sum
 
+    def test_lmap_via_nnc(self):
+        """Test correct lmaps generated indirectly via nnc_table.
+        Cells in nnc_table for layer number completely controlled by lmap.
+        """
+        grid, region, rid = _make_box_grid_with_region(dimension=(3, 3, 3))
+
+        region.values = np.ones(region.values.shape)
+        region.values[1][1][1] = rid
+
+        _, nnc_table = create_nested_hybrid_grid(
+            grid, region, rid, refinement=(2, 2, 2)
+        )
+        nnc_table1 = nnc_table[
+            (nnc_table["I1"] == 2)
+            & (nnc_table["J1"] == 1)
+            & (nnc_table["I2"] == 5)
+            & (nnc_table["J2"] == 1)
+        ]
+
+        assert nnc_table1[(nnc_table1["K1"] == 2) & (nnc_table1["K2"] == 2)].shape == (
+            1,
+            7,
+        )
+        assert nnc_table1[(nnc_table1["K1"] == 2) & (nnc_table1["K2"] == 3)].shape == (
+            1,
+            7,
+        )
+        assert nnc_table1[(nnc_table1["K1"] == 3) & (nnc_table1["K2"] == 2)].shape == (
+            0,
+            7,
+        )
+
+    def test_lmap_generation_simple(self):
+        """Tests that the correct layer mappings are generated"""
+
+        lmap1, lmap2 = _generate_layer_mappings(3, 2, (2, 2, 2), (1, 1, 1))
+
+        assert np.array_equal(lmap1, np.array([0, 1, 3]))
+        assert np.array_equal(lmap2, np.array([1, 2]))
+
+    def test_lmap_generation_no_offset(self):
+        """Tests that the correct layer mappings are generated"""
+
+        lmap1, lmap2 = _generate_layer_mappings(3, 4, (2, 2, 2), (1, 1, 0))
+
+        assert np.array_equal(lmap1, np.array([0, 2, 4]))
+        assert np.array_equal(lmap2, np.array([0, 1, 2, 3]))
+
+    def test_lmap_generation_full_offset(self):
+        """Tests that the correct layer mappings are generated"""
+
+        lmap1, lmap2 = _generate_layer_mappings(3, 4, (2, 2, 2), (1, 1, 3))
+
+        assert np.array_equal(lmap1, np.array([0, 1, 2]))
+        assert np.array_equal(lmap2, np.array([3, 4, 5, 6]))
+
+    def test_lmap_generation_ref10(self):
+        """Tests that the correct layer mappings are generated"""
+
+        lmap1, lmap2 = _generate_layer_mappings(3, 20, (2, 2, 10), (1, 1, 1))
+
+        assert np.array_equal(lmap1, np.array([0, 1, 11]))
+        assert np.array_equal(lmap2, np.arange(20) + 1)
+
 
 # ---------------------------------------------------------------------------
 # Tests for get_transmissibilities with nested hybrid NNCs
 # ---------------------------------------------------------------------------
 
 
 class TestTransmissibilitiesOnMergedGrid:
-    """Test calling get_transmissibilities on the merged grid with NEST_ID."""
+    """Test calling get_transmissibilities on the merged grid"""
 
     @staticmethod
     def _build_merged_with_props(dimension=(6, 6, 2), refinement=(1, 1, 1)):