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
@@ -52,8 +51,8 @@ class StructuredDiscretization(Package, IRegridPackage, IMaskingSettings):
         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
+        solution. If the idomain value for a cell is >0, the cell exists in the
+        simulation. If the idomain value for a cell is <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.
@@ -214,7 +213,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,19 @@ 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 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 >0, the cell exists in the simulation. If the idomain
+        value for a cell is <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.
     validate: {True, False}
         Flag to indicate whether the package should be validated upon
         initialization. This raises a ValidationError if package input is
@@ -57,7 +68,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,7 +325,7 @@ 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
+    * lower "pass-through" cells (idomain <0) are multitude of (number of
       cells in a layer)
 
     Parameters

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,12 +17,12 @@ 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
     idomain_float = idomain_float.combine_first(

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):