Issue #1376 active idomain larger than zero

docs/api/changelog.rst

@@ -23,6 +23,10 @@ Changed
 
 - :func:`imod.formats.prj.open_projectfile_data` now also assigns negative and
   zero layer numbers to grid coordinates.
+- In :class:`imod.mf6.StructuredDiscretization`, IDOMAIN can now respectively be
+  > 0 to indicate an active cell and <0 to indicate a vertical passthrough cell,
+  consistent with MODFLOW6. Previously this could only be indicated with 1 and
+  -1.
 
 
 [1.0.0rc1] - 2024-12-20

imod/mf6/clipped_boundary_condition_creator.py

@@ -28,7 +28,7 @@ def create_clipped_boundary(
         covered by other ConstantHead packages
 
     """
-    active_grid_boundary = active_grid_boundary_xy(idomain == 1)
+    active_grid_boundary = active_grid_boundary_xy(idomain > 0)
     unassigned_grid_boundaries = _find_unassigned_grid_boundaries(
         active_grid_boundary, original_constant_head_boundaries
     )

imod/mf6/dis.py

@@ -27,7 +27,6 @@
     DTypeSchema,
     IdentityNoDataSchema,
     IndexesSchema,
-    UniqueValuesSchema,
     ValidationError,
 )
 from imod.typing.grid import GridDataArray, is_unstructured
@@ -47,16 +46,20 @@ class StructuredDiscretization(Package, IRegridPackage, IMaskingSettings):
     bottom: array of floats (xr.DataArray)
         is the bottom elevation for each cell.
     idomain: array of integers (xr.DataArray)
-        Indicates the existence status of a cell. Horizontal discretization
-        information will be derived from the x and y coordinates of the
-        DataArray. If the idomain value for a cell is 0, the cell does not exist
-        in the simulation. Input and output values will be read and written for
-        the cell, but internal to the program, the cell is excluded from the
-        solution. If the idomain value for a cell is 1, the cell exists in the
-        simulation. if the idomain value for a cell is -1, the cell does not
-        exist in the simulation. Furthermore, the first existing cell above will
-        be connected to the first existing cell below. This type of cell is
-        referred to as a "vertical pass through" cell.
+        Indicates the existence status of a cell as follows:
+
+        * If 0, the cell does not exist in the simulation. Input and output
+          values will be read and written for the cell, but internal to the
+          program, the cell is excluded from the solution.
+        * If >0, the cell exists in the simulation.
+        * If <0, the cell does not exist in the simulation. Furthermore, the
+          first existing cell above will be connected to the first existing cell
+          below. This type of cell is referred to as a "vertical pass through"
+          cell.
+
+        This DataArray needs contain a ``"layer"``, ``"y"`` and ``"x"``
+        coordinate. Horizontal discretization information will be derived from
+        its x and y coordinates.
     validate: {True, False}
         Flag to indicate whether the package should be validated upon
         initialization. This raises a ValidationError if package input is
@@ -214,7 +217,6 @@ def from_imod5_data(
 
         # Validate iMOD5 data
         if validate:
-            UniqueValuesSchema([-1, 0, 1]).validate(imod5_data["bnd"]["ibound"])
             if not np.all(
                 new_package_data["top"][1:].data == new_package_data["bottom"][:-1].data
             ):

imod/mf6/disv.py

@@ -27,8 +27,24 @@ class VerticesDiscretization(Package, IRegridPackage, IMaskingSettings):
     Parameters
     ----------
     top: array of floats (xu.UgridDataArray)
+        is the top elevation for each cell in the top model layer.
     bottom: array of floats (xu.UgridDataArray)
+        is the bottom elevation for each cell.
     idomain: array of integers (xu.UgridDataArray)
+        Indicates the existence status of a cell.
+
+        * If 0, the cell does not exist in the simulation. Input and output
+          values will be read and written for the cell, but internal to the
+          program, the cell is excluded from the solution.
+        * If >0, the cell exists in the simulation.
+        * If <0, the cell does not exist in the simulation. Furthermore, the
+          first existing cell above will be connected to the first existing cell
+          below. This type of cell is referred to as a "vertical pass through"
+          cell.
+
+        This UgridDataArray needs to contain a ``"layer"`` coordinate and a face
+        dimension. Horizontal discretization information will be derived from
+        its face dimension.
     validate: {True, False}
         Flag to indicate whether the package should be validated upon
         initialization. This raises a ValidationError if package input is
@@ -57,7 +73,7 @@ class VerticesDiscretization(Package, IRegridPackage, IMaskingSettings):
     _write_schemata = {
         "idomain": (AnyValueSchema(">", 0),),
         "top": (
-            AllValueSchema(">", "bottom", ignore=("idomain", "==", -1)),
+            AllValueSchema(">", "bottom", ignore=("idomain", "<=", 0)),
             IdentityNoDataSchema(other="idomain", is_other_notnull=(">", 0)),
             # No need to check coords: dataset ensures they align with idomain.
         ),

imod/mf6/hfb.py

@@ -557,7 +557,8 @@ def _to_connected_cells_dataset(
             - value name
         """
         top, bottom = broadcast_to_full_domain(idomain, top, bottom)
-        k = idomain * k
+        # Broadcast to grid (if necessary) and deactivate inactive cells.
+        k = (idomain > 0) * k
         # Enforce unstructured
         unstructured_grid, top, bottom, k = (
             as_ugrid_dataarray(grid) for grid in [idomain, top, bottom, k]

imod/mf6/model.py

@@ -672,18 +672,18 @@ def mask_all_packages(
         """
         This function applies a mask to all packages in a model. The mask must
         be presented as an idomain-like integer array that has 0 (inactive) or
-        -1 (vertical passthrough) values in filtered cells and 1 in active
+        <0 (vertical passthrough) values in filtered cells and >0 in active
         cells.
-        Masking will overwrite idomain with the mask where the mask is 0 or -1.
-        Where the mask is 1, the original value of idomain will be kept. Masking
+        Masking will overwrite idomain with the mask where the mask is <=0.
+        Where the mask is >0, the original value of idomain will be kept. Masking
         will update the packages accordingly, blanking their input where needed,
         and is therefore not a reversible operation.
 
         Parameters
         ----------
         mask: xr.DataArray, xu.UgridDataArray of ints
-            idomain-like integer array. 1 sets cells to active, 0 sets cells to inactive,
-            -1 sets cells to vertical passthrough
+            idomain-like integer array. >0 sets cells to active, 0 sets cells to inactive,
+            <0 sets cells to vertical passthrough
         """
 
         _mask_all_packages(self, mask)

imod/mf6/multimodel/exchange_creator_unstructured.py

@@ -104,7 +104,7 @@ def _create_global_to_local_idx(
             local_cell_indices = cls._get_local_cell_indices(submodel_partition_info)
 
             global_cell_indices_partition = global_cell_indices.where(
-                submodel_partition_info.active_domain == 1
+                submodel_partition_info.active_domain > 0
             )
             global_cell_indices_partition = global_cell_indices_partition.dropna(
                 "mesh2d_nFaces", how="all"

imod/mf6/out/dis.py

@@ -325,8 +325,8 @@ def dis_indices(
     * the right neighbor is +(1)
     * the front neighbor is +(number of cells in a column)
     * the lower neighbor is +(number of cells in a layer)
-    * lower "pass-through" cells (idomain == -1) are multitude of (number of
-      cells in a layer)
+    * lower "vertical passthrough" cells (idomain <0) are multitude of (number
+      of cells in a layer)
 
     Parameters
     ----------

imod/mf6/out/disv.py

@@ -436,8 +436,8 @@ def disv_indices(
 
     # Vertical flows
     # --------------
-    # Stored in an array of shape (nlayer, nface).
-    # For pass-through cells, set it from layers i to j using ragged_arange.
+    # Stored in an array of shape (nlayer, nface). For vertical passthrough
+    # cells, set it from layers i to j using ragged_arange.
     n_pass = diff[is_vertical] // ncells_per_layer
     ii = np.repeat(i[is_vertical], n_pass) + ragged_arange(n_pass) * ncells_per_layer
     lower.ravel()[ii] = np.repeat(index[is_vertical], n_pass)

imod/mf6/package.py

@@ -565,8 +565,8 @@ def mask(self, mask: GridDataArray) -> Any:
         Parameters
         ----------
         mask: xr.DataArray, xu.UgridDataArray of ints
-            idomain-like integer array. 1 sets cells to active, 0 sets cells to inactive,
-            -1 sets cells to vertical passthrough
+            idomain-like integer array. >0 sets cells to active, 0 sets cells to inactive,
+            <0 sets cells to vertical passthrough
 
         Returns
         -------

imod/mf6/simulation.py

@@ -1335,16 +1335,16 @@ def mask_all_models(
         """
         This function applies a mask to all models in a simulation, provided they use
         the same discretization. The  method parameter "mask" is an idomain-like array.
-        Masking will overwrite idomain with the mask where the mask is 0 or -1.
-        Where the mask is 1, the original value of idomain will be kept.
+        Masking will overwrite idomain with the mask where the mask is <0.
+        Where the mask is >0, the original value of idomain will be kept.
         Masking will update the packages accordingly, blanking their input where needed,
         and is therefore not a reversible operation.
 
         Parameters
         ----------
         mask: xr.DataArray, xu.UgridDataArray of ints
-            idomain-like integer array. 1 sets cells to active, 0 sets cells to inactive,
-            -1 sets cells to vertical passthrough
+            idomain-like integer array. >0 sets cells to active, 0 sets cells to inactive,
+            <0 sets cells to vertical passthrough
         """
         _mask_all_models(self, mask)
 

imod/mf6/utilities/grid.py

@@ -39,9 +39,10 @@ def broadcast_to_full_domain(
     """
     Broadcast the bottom and top array to have the same shape as the idomain
     """
-    bottom = idomain * bottom
+    active = idomain > 0
+    bottom = active * bottom
     top = (
-        idomain * top
+        active * top
         if hasattr(top, "coords") and "layer" in top.coords
         else create_layered_top(bottom, top)
     )

imod/mf6/utilities/imod5_converter.py

@@ -17,14 +17,14 @@ def convert_ibound_to_idomain(
     idomain = np.abs(ibound)
 
     # Thickness <= 0 -> IDOMAIN = -1
-    active_and_zero_thickness = (thickness <= 0) & (idomain == 1)
+    active_and_zero_thickness = (thickness <= 0) & (idomain > 0)
     # Don't make cells at top or bottom vpt, these should be inactive.
     # First, set all potential vpts to nan to be able to utilize ffill and bfill
     idomain_float = idomain.where(~active_and_zero_thickness)  # type: ignore[attr-defined]
-    passthrough = (idomain_float.ffill("layer") == 1) & (
-        idomain_float.bfill("layer") == 1
+    passthrough = (idomain_float.ffill("layer") > 0) & (
+        idomain_float.bfill("layer") > 0
     )
-    # Then fill nans where passthrough with -1
+    # Then fill nans where vertical passthrough with -1
     idomain_float = idomain_float.combine_first(
         full_like(idomain_float, -1.0, dtype=float).where(passthrough)
     )

imod/mf6/utilities/mask.py

@@ -99,14 +99,15 @@ def _skip_variable(package: IMaskingSettings, var: str) -> bool:
 def _mask_spatial_var(self, var: str, mask: GridDataArray) -> GridDataArray:
     da = self.dataset[var]
     array_mask = _adjust_mask_for_unlayered_data(da, mask)
+    active = array_mask > 0
 
     if issubclass(da.dtype.type, numbers.Integral):
         if var == "idomain":
-            return da.where(array_mask > 0, other=array_mask)
+            return da.where(active, other=array_mask)
         else:
-            return da.where(array_mask > 0, other=0)
+            return da.where(active, other=0)
     elif issubclass(da.dtype.type, numbers.Real):
-        return da.where(array_mask > 0)
+        return da.where(active)
     else:
         raise TypeError(
             f"Expected dtype float or integer. Received instead: {da.dtype}"

imod/mf6/wel.py

@@ -73,7 +73,7 @@ def _assign_dims(arg: Any) -> Tuple | xr.DataArray:
 def mask_2D(package: GridAgnosticWell, domain_2d: GridDataArray) -> GridAgnosticWell:
     point_active = points_values(domain_2d, x=package.x, y=package.y)
 
-    is_inside_exterior = point_active == 1
+    is_inside_exterior = point_active > 0
     selection = package.dataset.loc[{"index": is_inside_exterior}]
 
     cls = type(package)

imod/msw/utilities/imod5_converter.py

@@ -83,9 +83,7 @@ def is_msw_active_cell(
         Cells active per subunit
     """
     mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True)
-    subunit_active = (
-        (imod5_cap["boundary"] == 1) & (msw_area > 0) & (mf6_top_active >= 1)
-    )
+    subunit_active = (imod5_cap["boundary"] > 0) & (msw_area > 0) & (mf6_top_active > 0)
     active = subunit_active.any(dim="subunit")
     return MetaSwapActive(active, subunit_active)
 

imod/prepare/cleanup.py

@@ -41,7 +41,7 @@ def _cleanup_robin_boundary(
     idomain: GridDataArray, grids: dict[str, GridDataArray]
 ) -> dict[str, GridDataArray]:
     """Cleanup robin boundary condition (i.e. bc with conductance)"""
-    active = idomain == 1
+    active = idomain > 0
     # Deactivate conductance cells outside active domain; this nodata
     # inconsistency will be aligned in the final call to align_nodata
     conductance = grids["conductance"].where(active)

imod/prepare/topsystem/allocation.py

@@ -83,7 +83,7 @@ def allocate_riv_cells(
         enumerator.
     active: DataArray | UgridDatarray
         Boolean array containing active model cells. For Modflow 6, this is the
-        equivalent of ``idomain == 1``.
+        equivalent of ``idomain > 0``.
     top: DataArray | UgridDatarray
         Grid containing tops of model layers. If has no layer dimension, is
         assumed as top of upper layer and the other layers are padded with
@@ -156,7 +156,7 @@ def allocate_drn_cells(
         enumerator.
     active: DataArray | UgridDatarray
         Boolean array containing active model cells. For Modflow 6, this is the
-        equivalent of ``idomain == 1``.
+        equivalent of ``idomain > 0``.
     top: DataArray | UgridDatarray
         Grid containing tops of model layers. If has no layer dimension, is
         assumed as top of upper layer and the other layers are padded with
@@ -217,7 +217,7 @@ def allocate_ghb_cells(
         enumerator.
     active: DataArray | UgridDatarray
         Boolean array containing active model cells. For Modflow 6, this is the
-        equivalent of ``idomain == 1``.
+        equivalent of ``idomain > 0``.
     top: DataArray | UgridDatarray
         Grid containing tops of model layers. If has no layer dimension, is
         assumed as top of upper layer and the other layers are padded with
@@ -275,7 +275,7 @@ def allocate_rch_cells(
         enumerator.
     active: DataArray | UgridDataArray
         Boolean array containing active model cells. For Modflow 6, this is the
-        equivalent of ``idomain == 1``.
+        equivalent of ``idomain > 0``.
     rate: DataArray | UgridDataArray
         Array with recharge rates. This will only be used to infer where
         recharge cells are defined.

imod/tests/fixtures/backward_compatibility_fixture.py

@@ -21,7 +21,7 @@ def imod5_dataset():
 
     # Fix data for ibound  as it contains floating values like 0.34, 0.25 etc.
     ibound = data["bnd"]["ibound"]
-    ibound = ibound.where(ibound <= 0, 1)
+    ibound = ibound.where(ibound <= 0, 2)
     data["bnd"]["ibound"] = ibound
     return data, pd
 

imod/tests/fixtures/mf6_small_models_fixture.py

@@ -123,7 +123,7 @@ def structured_flow_model() -> imod.mf6.GroundwaterFlowModel:
     gwf_model = _make_model(grid_data_structured, cellsize)
 
     bottom = grid_data_structured_layered(np.float64, -1.0, cellsize)
-    idomain = grid_data_structured(np.int32, 1, cellsize)
+    idomain = grid_data_structured(np.int32, 2, cellsize)
 
     # Unassign the constant-head package from some cells to have a more meaningful simulation
     constant_head = gwf_model["chd"].dataset["head"]
@@ -152,7 +152,7 @@ def unstructured_flow_model() -> imod.mf6.GroundwaterFlowModel:
     constant_head.loc[{"mesh2d_nFaces": 13, "layer": 2}] = np.nan
 
     bottom = grid_data_unstructured_layered(np.float64, -1.0, cellsize)
-    idomain = grid_data_unstructured(np.int32, 1, cellsize)
+    idomain = grid_data_unstructured(np.int32, 2, cellsize)
     gwf_model["disv"] = imod.mf6.VerticesDiscretization(
         top=10.0, bottom=bottom, idomain=idomain
     )

imod/tests/test_mf6/test_mf6_dis.py

@@ -209,7 +209,9 @@ def test_from_imod5_data__idomain_values(imod5_dataset):
     # Test if idomain has appropriate count
     assert (dis["idomain"] == -1).sum() == 371824
     assert (dis["idomain"] == 0).sum() == 176912
-    assert (dis["idomain"] == 1).sum() == 703936
+    # IBOUND -1 was converted to 1, therefore not all active cells are 2
+    assert (dis["idomain"] == 1).sum() == 15607
+    assert (dis["idomain"] == 2).sum() == 688329
 
 
 @pytest.mark.usefixtures("imod5_dataset")
@@ -256,6 +258,8 @@ def test_from_imod5_data__validation_error(tmp_path):
     tmp_path = imod.util.temporary_directory()
     data = imod.data.imod5_projectfile_data(tmp_path)
     data = data[0]
+    # Set bottom above top, to create "gap" intebetween interfaces.
+    data["bot"]["bottom"][20, 6, 6] = data["top"]["top"][21, 6, 6] - 1.0
 
     _load_imod5_data_in_memory(data)
     with pytest.raises(ValidationError):

imod/tests/test_mf6/test_mf6_model.py

@@ -432,4 +432,4 @@ def test_get_domain_geometry(structured_flow_model: GroundwaterFlowModel):
 
     assert np.all(np.isin(top.values, [10.0]))
     assert np.all(np.isin(bottom.values, [-1.0, -2.0, -3.0]))
-    assert np.all(np.isin(idomain.values, [1]))
+    assert np.all(np.isin(idomain.values, [2]))

imod/tests/test_mf6/test_mf6_regrid_model.py

@@ -100,7 +100,7 @@ def test_regrid_unstructured_model_with_inactive_cells(
     assert not validation_result.has_errors()
     new_idomain = new_gwf_model.domain
     assert (
-        new_idomain.max().values[()] == 1
+        new_idomain.max().values[()] == 2
         and new_idomain.min().values[()] == inactivity_marker
     )
 

imod/tests/test_mf6/test_utilities/test_imod5_converter.py

@@ -62,6 +62,12 @@ def case_mixed(self):
         idomain = [1, -1, 1, 0, 0]
         return thickness, ibound, idomain
 
+    def case_mixed__idomain_2(self):
+        thickness = [1.0, 0.0, 1.0, 0.0, 1.0]
+        ibound = [2, -1, 2, 2, 0]
+        idomain = [2, -1, 2, 0, 0]
+        return thickness, ibound, idomain
+
 
 @parametrize_with_cases(argnames="thickness,ibound,expected", cases=IboundCases)
 def test_convert_ibound_to_idomain(template_grid, thickness, ibound, expected):