From c3558b76a4dd2dd75fc1a47d9bf81ff712ff754b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 22 Sep 2026 14:30:48 +0200 Subject: [PATCH 01/39] Add mask_packages method --- docs/api/changelog.rst | 4 +++ docs/api/mf6.rst | 2 ++ imod/common/interfaces/imodel.py | 14 +++++++++- imod/common/utilities/mask.py | 13 +++++---- imod/mf6/model.py | 46 ++++++++++++++++++++++++++++---- 5 files changed, 68 insertions(+), 11 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index a5b40a49a..6f02e203c 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -18,6 +18,10 @@ Added :meth:`imod.msw.SprinklingPoints.from_imod5_data`. - :class:`imod.mf6.LayeredWell.from_imod5_cap_data` now also supports loading wells from IPF files in an iMOD5 CAP dataset. +- Added :meth:`imod.mf6.GroundwaterFlowModel.mask_packages` and + :meth:`imod.mf6.GroundwaterTransportModel.mask_packages` to mask specific + packages of a groundwater flow model and a groundwater transport model + respectively. Fixed ~~~~~ diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst index 32fb67246..dda794721 100644 --- a/docs/api/mf6.rst +++ b/docs/api/mf6.rst @@ -46,6 +46,7 @@ Model objects & methods Modflow6Simulation.set_validation_settings GroundwaterFlowModel GroundwaterFlowModel.mask_all_packages + GroundwaterFlowModel.mask_packages GroundwaterFlowModel.prepare_wel_for_mf6 GroundwaterFlowModel.regrid_like GroundwaterFlowModel.dump @@ -60,6 +61,7 @@ Model objects & methods GroundwaterFlowModel.get_diskey GroundwaterTransportModel GroundwaterTransportModel.mask_all_packages + GroundwaterTransportModel.mask_packages GroundwaterTransportModel.dump GroundwaterTransportModel.clip_box GroundwaterTransportModel.regrid_like diff --git a/imod/common/interfaces/imodel.py b/imod/common/interfaces/imodel.py index 2fe91e0a4..0b30e0a6e 100644 --- a/imod/common/interfaces/imodel.py +++ b/imod/common/interfaces/imodel.py @@ -16,9 +16,21 @@ class IModel(IDict): def mask_all_packages(self, mask: GridDataArray): raise NotImplementedError + @abstractmethod + def mask_packages( + self, + package_names: list[str], + mask: GridDataArray, + ignore_time_purge_empty: bool = False, + ): + raise NotImplementedError + @abstractmethod def purge_empty_packages( - self, model_name: Optional[str] = "", ignore_time: bool = False + self, + model_name: Optional[str] = "", + ignore_time: bool = False, + package_names: list[str] | None = None, ) -> None: raise NotImplementedError diff --git a/imod/common/utilities/mask.py b/imod/common/utilities/mask.py index 21f689f27..eba69d86f 100644 --- a/imod/common/utilities/mask.py +++ b/imod/common/utilities/mask.py @@ -62,15 +62,18 @@ def mask_all_models( ) -def mask_all_packages( +def mask_packages( model: IModel, + package_names: list[str], mask: GridDataArray, ignore_time_purge_empty: bool = False, -): +) -> None: _validate_coords_mask(mask) - for pkgname, pkg in model.items(): - model[pkgname] = pkg.mask(mask) - model.purge_empty_packages(ignore_time=ignore_time_purge_empty) + for pkgname in package_names: + model[pkgname] = model[pkgname].mask(mask) + model.purge_empty_packages( + ignore_time=ignore_time_purge_empty, package_names=package_names + ) def mask_package(package: IPackage, mask: GridDataArray) -> IPackage: diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 2a931732f..54208cefe 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -22,7 +22,7 @@ from imod.common.statusinfo import NestedStatusInfo, StatusInfo, StatusInfoBase from imod.common.utilities.clip import clip_box_dataset from imod.common.utilities.dump_model import dump_model -from imod.common.utilities.mask import mask_all_packages +from imod.common.utilities.mask import mask_packages from imod.common.utilities.regrid import _regrid_like from imod.common.utilities.schemata import ( concatenate_schemata_dicts, @@ -930,11 +930,41 @@ def mask_all_packages( Whether to ignore time dimension when purging empty packages. Can improve performance when masking models with many time steps. """ + package_names = list(self.keys()) + mask_packages(self, package_names, mask, ignore_time_purge_empty) - mask_all_packages(self, mask, ignore_time_purge_empty) + def mask_packages( + self, + package_names: list[str], + mask: GridDataArray, + ignore_time_purge_empty: bool = False, + ) -> None: + """ + This function applies a mask to packages in a model. The mask must + be presented as an idomain-like integer array that has 0 (inactive) or + <0 (vertical passthrough) values in filtered cells and >0 in active + cells. + 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. >0 sets cells to active, 0 sets cells to inactive, + <0 sets cells to vertical passthrough + ignore_time_purge_empty: bool, default False + Whether to ignore time dimension when purging empty packages. Can + improve performance when masking models with many time steps. + """ + mask_packages(self, package_names, mask, ignore_time_purge_empty) def purge_empty_packages( - self, model_name: Optional[str] = "", ignore_time: bool = False + self, + model_name: Optional[str] = "", + ignore_time: bool = False, + package_names: Optional[list[str]] = None, ) -> None: """ This method removes empty packages from the model in place. @@ -948,11 +978,17 @@ def purge_empty_packages( timesteps. If True, packages are considered empty if they have no data at the first time step. The latter can increase performance considerably. + package_names: list[str], optional + List of package names to check for emptiness. If None, all packages + are checked. """ + if package_names is None: + package_names = list(self.keys()) + empty_packages = [ package_name - for package_name, package in self.items() - if package.is_empty(ignore_time=ignore_time) + for package_name in package_names + if self[package_name].is_empty(ignore_time=ignore_time) ] logger.info( f"packages: {empty_packages} removed in {model_name}, because all empty" From 1c8ffb0e6d3ac9f0f6bc141fa4cb7aeddf0fe379 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 22 Sep 2026 15:30:32 +0200 Subject: [PATCH 02/39] Add mask_topsystem_packages utility function --- imod/mf6/model_gwf.py | 8 ++++++ imod/mf6/utilities/imod5_converter.py | 37 ++++++++++++++++++++++++++- 2 files changed, 44 insertions(+), 1 deletion(-) diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 545f948da..7e5d764af 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -37,6 +37,7 @@ from imod.mf6.riv import River from imod.mf6.sto import StorageCoefficient from imod.mf6.utilities.chd_concat import concat_layered_chd_packages +from imod.mf6.utilities.imod5_converter import mask_topsystem_packages from imod.mf6.validation_settings import ValidationSettings from imod.mf6.wel import LayeredWell, Well from imod.prepare.topsystem.default_allocation_methods import ( @@ -453,4 +454,11 @@ def from_imod5_data( for key, chd_package in chd_packages.items(): result[key] = chd_package + # Mask all topsystem packages where IBOUND == -1 + mask_topsystem_packages( + imod5_data, + result, + cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")), + regrid_cache, + ) return result diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index e27c21581..eb3a19bb6 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -4,10 +4,12 @@ import pandas as pd import xarray as xr +from imod.common.interfaces.imodel import IModel from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data from imod.mf6.package import Package +from imod.mf6.regrid.regrid_schemes import ConstantHeadRegridMethod from imod.typing import GridDataArray, GridDataDict, Imod5DataDict from imod.typing.grid import full_like from imod.util.regrid import RegridderWeightsCache @@ -140,7 +142,7 @@ def well_from_imod5_cap_data( def regrid_imod5_pkg_data( - cls: type[Package], + cls: Optional[type[Package]], imod5_pkg_data: GridDataDict, target_dis: Package, regridder_types: Optional[DataclassType], @@ -150,6 +152,11 @@ def regrid_imod5_pkg_data( Regrid iMOD5 package data to target idomain. Optionally get regrid methods from class if not provided. """ + if cls is None and regridder_types is None: + raise ValueError( + "Either cls or regridder_types must be provided for regridding." + ) + target_idomain = target_dis.dataset["idomain"] # set up regridder methods @@ -175,3 +182,31 @@ def chd_cells_from_imod5_data( head = head.where(target_idomain > 0) return {"head": head} + + +def mask_topsystem_packages( + imod5_data: Imod5DataDict, + model: IModel, + regridder_types: ConstantHeadRegridMethod, + regrid_cache: RegridderWeightsCache, +) -> None: + """ + Mask all top system packages where IBOUND == -1. + """ + from imod.mf6.topsystem import TopSystemBoundaryCondition + + ibound = imod5_data["bnd"]["ibound"] + regridded_ibound = regrid_imod5_pkg_data( + cls=None, + imod5_pkg_data={"ibound": ibound}, + target_dis=model["dis"], + regridder_types=regridder_types, + regrid_cache=regrid_cache, + )["ibound"] + mask = regridded_ibound == -1 + + topsystem_packages = [ + key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition) + ] + for key in topsystem_packages: + model[key].mask(mask) From be94e0f32a30f88d1daa18aded8fcef6a5b137bb Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 12:36:15 +0200 Subject: [PATCH 03/39] Provide proper mask and fix mypy issues --- imod/mf6/utilities/imod5_converter.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index eb3a19bb6..31995448d 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -1,4 +1,4 @@ -from typing import Optional, Union +from typing import Optional, Union, cast import numpy as np import pandas as pd @@ -152,16 +152,18 @@ def regrid_imod5_pkg_data( Regrid iMOD5 package data to target idomain. Optionally get regrid methods from class if not provided. """ - if cls is None and regridder_types is None: + if (cls is None) and (regridder_types is None): raise ValueError( "Either cls or regridder_types must be provided for regridding." ) + # set up regridder methods + elif (cls is not None) and (regridder_types is None): # check cls not None for mypy + regridder_types = cls.get_regrid_methods() + # For mypy to succeed + regridder_types = cast(DataclassType, regridder_types) target_idomain = target_dis.dataset["idomain"] - # set up regridder methods - if regridder_types is None: - regridder_types = cls.get_regrid_methods() # regrid the input data regridded_pkg_data = _regrid_package_data( imod5_pkg_data, target_idomain, regridder_types, regrid_cache, {} @@ -185,16 +187,21 @@ def chd_cells_from_imod5_data( def mask_topsystem_packages( - imod5_data: Imod5DataDict, + imod5_data: dict[str, dict[str, GridDataArray]], model: IModel, - regridder_types: ConstantHeadRegridMethod, + regridder_types: Optional[ConstantHeadRegridMethod], regrid_cache: RegridderWeightsCache, ) -> None: """ - Mask all top system packages where IBOUND == -1. + Mask all top system packages where IBOUND < 0. These locations are assigned + a constant head. """ + # Import here to avoid circular import issues from imod.mf6.topsystem import TopSystemBoundaryCondition + if regridder_types is None: + regridder_types = ConstantHeadRegridMethod() + ibound = imod5_data["bnd"]["ibound"] regridded_ibound = regrid_imod5_pkg_data( cls=None, @@ -203,10 +210,10 @@ def mask_topsystem_packages( regridder_types=regridder_types, regrid_cache=regrid_cache, )["ibound"] - mask = regridded_ibound == -1 + is_active = regridded_ibound >= 0 topsystem_packages = [ key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition) ] for key in topsystem_packages: - model[key].mask(mask) + model[key] = model[key].mask(is_active) From 69bd29c852791bad0e35ab188afefdec32759991 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 12:37:35 +0200 Subject: [PATCH 04/39] Add test --- imod/tests/test_mf6/test_mf6_simulation.py | 37 ++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py index e980b0eb9..3c1c1a718 100644 --- a/imod/tests/test_mf6/test_mf6_simulation.py +++ b/imod/tests/test_mf6/test_mf6_simulation.py @@ -610,6 +610,43 @@ def test_import_from_imod5(imod5_dataset, tmp_path): assert simulation._validation_context.strict_well_validation is False +@pytest.mark.unittest_jit +def test_import_from_imod5_mask_topsystem(imod5_dataset): + """Test importing from imod5 masks the top system packages""" + # Arrange + imod5_data = imod5_dataset[0] + period_data = imod5_dataset[1] + + datelist = pd.date_range(start="1/1/1989", end="1/1/2013", freq="W") + # Act + simulation = Modflow6Simulation.from_imod5_data( + imod5_data, + period_data, + datelist, + SimulationAllocationOptions, + SimulationDistributingOptions, + ) + # Assert + ibound = imod5_data["bnd"]["ibound"].isel(layer=0, drop=True) + topsystem_mask = ibound < 0 + topsystem_keys = [ + "rch", + "drn-1", + "drn-2", + "riv-1riv", + "riv-1drn", + "riv-2riv", + "riv-2drn", + ] + for key in topsystem_keys: + topsystem_pkg = simulation["imported_model"][key] + # Take first grid var + gridded_var = topsystem_pkg.dataset[topsystem_pkg._period_data[0]] + # True wherever topsystem is masked but gridded_var still has a value + bad = gridded_var.notnull() & topsystem_mask + assert not bad.any().item() + + @pytest.mark.unittest_jit def test_import_from_imod5__custom_name(imod5_dataset): imod5_data = imod5_dataset[0] From d0970c3d0facba5bf1873840237472c359351ef3 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 12:37:42 +0200 Subject: [PATCH 05/39] Fix docstring --- imod/mf6/model_gwf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 7e5d764af..0faf210eb 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -454,7 +454,7 @@ def from_imod5_data( for key, chd_package in chd_packages.items(): result[key] = chd_package - # Mask all topsystem packages where IBOUND == -1 + # Mask all topsystem packages where IBOUND < 0 mask_topsystem_packages( imod5_data, result, From c8cbc5f2b029531da0f874518e2fe802f29c9355 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 13:04:11 +0200 Subject: [PATCH 06/39] Move mask topsystem function to separate utility and add ITopSystemBoundaryCondition interface. --- imod/common/interfaces/itopsystembc.py | 19 +++++++++++++++++++ imod/mf6/model_gwf.py | 4 ++-- imod/mf6/topsystem.py | 5 ++++- imod/mf6/utilities/imod5_converter.py | 11 +++-------- imod/mf6/utilities/mask.py | 25 +++++++++++++++++++++++++ 5 files changed, 53 insertions(+), 11 deletions(-) create mode 100644 imod/common/interfaces/itopsystembc.py create mode 100644 imod/mf6/utilities/mask.py diff --git a/imod/common/interfaces/itopsystembc.py b/imod/common/interfaces/itopsystembc.py new file mode 100644 index 000000000..e202064c4 --- /dev/null +++ b/imod/common/interfaces/itopsystembc.py @@ -0,0 +1,19 @@ +from abc import abstractmethod + +from imod.common.interfaces.ipackage import IPackage +from imod.typing import GridDataDict, GridDataset + + +class ITopSystemBoundaryCondition(IPackage): + """ + Interface for top system boundary condition packages in MODFLOW 6. + """ + + @classmethod + @abstractmethod + def aggregate_layers(cls, dataset: GridDataset) -> GridDataDict: + raise NotImplementedError + + @abstractmethod + def reallocate(self): + raise NotImplementedError diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 0faf210eb..5296bdc62 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -37,7 +37,7 @@ from imod.mf6.riv import River from imod.mf6.sto import StorageCoefficient from imod.mf6.utilities.chd_concat import concat_layered_chd_packages -from imod.mf6.utilities.imod5_converter import mask_topsystem_packages +from imod.mf6.utilities.imod5_converter import mask_topsystem_packages_with_ibound from imod.mf6.validation_settings import ValidationSettings from imod.mf6.wel import LayeredWell, Well from imod.prepare.topsystem.default_allocation_methods import ( @@ -455,7 +455,7 @@ def from_imod5_data( result[key] = chd_package # Mask all topsystem packages where IBOUND < 0 - mask_topsystem_packages( + mask_topsystem_packages_with_ibound( imod5_data, result, cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")), diff --git a/imod/mf6/topsystem.py b/imod/mf6/topsystem.py index 14d67e0f3..692fa1d7a 100644 --- a/imod/mf6/topsystem.py +++ b/imod/mf6/topsystem.py @@ -3,6 +3,7 @@ from dataclasses import asdict from typing import Optional, Self, cast +from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition from imod.common.utilities.dataclass_type import DataclassType from imod.mf6.aggregate.aggregate_schemes import EmptyAggregationMethod from imod.mf6.boundary_condition import BoundaryCondition @@ -41,7 +42,9 @@ def _handle_reallocate_arguments( return allocation_option, distributing_option -class TopSystemBoundaryCondition(BoundaryCondition, abc.ABC): +class TopSystemBoundaryCondition( + BoundaryCondition, ITopSystemBoundaryCondition, abc.ABC +): """ Base class to add some extra functionality for topsystem packages, such as RCH, DRN, RIV, and GHB. diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index 31995448d..0da531a2a 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -10,6 +10,7 @@ from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data from imod.mf6.package import Package from imod.mf6.regrid.regrid_schemes import ConstantHeadRegridMethod +from imod.mf6.utilities.mask import mask_topsystem from imod.typing import GridDataArray, GridDataDict, Imod5DataDict from imod.typing.grid import full_like from imod.util.regrid import RegridderWeightsCache @@ -186,7 +187,7 @@ def chd_cells_from_imod5_data( return {"head": head} -def mask_topsystem_packages( +def mask_topsystem_packages_with_ibound( imod5_data: dict[str, dict[str, GridDataArray]], model: IModel, regridder_types: Optional[ConstantHeadRegridMethod], @@ -196,8 +197,6 @@ def mask_topsystem_packages( Mask all top system packages where IBOUND < 0. These locations are assigned a constant head. """ - # Import here to avoid circular import issues - from imod.mf6.topsystem import TopSystemBoundaryCondition if regridder_types is None: regridder_types = ConstantHeadRegridMethod() @@ -212,8 +211,4 @@ def mask_topsystem_packages( )["ibound"] is_active = regridded_ibound >= 0 - topsystem_packages = [ - key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition) - ] - for key in topsystem_packages: - model[key] = model[key].mask(is_active) + mask_topsystem(model, is_active) diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py new file mode 100644 index 000000000..574422d48 --- /dev/null +++ b/imod/mf6/utilities/mask.py @@ -0,0 +1,25 @@ +from imod.common.interfaces.imodel import IModel +from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition +from imod.typing import GridDataArray + + +def mask_topsystem(model: IModel, is_active: GridDataArray): + """ + Mask all top system packages in the model inplace with a boolean mask + indicating active cells. + + Parameters + ---------- + model : IModel + The MODFLOW 6 model containing top system packages. + is_active : GridDataArray + A boolean array indicating active cells. Top system packages will be masked + where this array is False. + """ + topsystem_packages = [ + key + for key, pkg in model.items() + if isinstance(pkg, ITopSystemBoundaryCondition) + ] + for key in topsystem_packages: + model[key] = model[key].mask(is_active) From 29204bd381952fd4a8ecf89282453f1bc4961980 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 13:09:28 +0200 Subject: [PATCH 07/39] Remove method --- imod/common/interfaces/itopsystembc.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/imod/common/interfaces/itopsystembc.py b/imod/common/interfaces/itopsystembc.py index e202064c4..de541ec09 100644 --- a/imod/common/interfaces/itopsystembc.py +++ b/imod/common/interfaces/itopsystembc.py @@ -13,7 +13,3 @@ class ITopSystemBoundaryCondition(IPackage): @abstractmethod def aggregate_layers(cls, dataset: GridDataset) -> GridDataDict: raise NotImplementedError - - @abstractmethod - def reallocate(self): - raise NotImplementedError From a4ba60d0ba03cfcc3f3deb338d8d2cc02f5d6970 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 14:02:25 +0200 Subject: [PATCH 08/39] Also clip topsystems for clip_box when states_for_boundary are provided --- imod/mf6/model.py | 8 ++++++++ imod/tests/test_mf6/test_ex01_twri.py | 11 +++++++++++ 2 files changed, 19 insertions(+) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 54208cefe..70f21325f 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -43,11 +43,13 @@ StateType, create_clipped_boundary, ) +from imod.mf6.utilities.mask import mask_topsystem from imod.mf6.utilities.mf6hfb import merge_hfb_packages from imod.mf6.validation_settings import ValidationSettings from imod.mf6.wel import GridAgnosticWell from imod.mf6.write_context import WriteContext from imod.schemata import SchemataDict, ValidationError +from imod.select.grid import active_grid_boundary_xy from imod.typing import GridDataArray from imod.util.regrid import RegridderWeightsCache @@ -813,6 +815,12 @@ def clip_box( if clipped_boundary_condition is not None: clipped[pkg_name] = clipped_boundary_condition + # Clip topsystem packages where the active grid boundary has + # changed. + _, _, idomain_clipped = clipped._get_domain_geometry() + active_bounds_clipped = active_grid_boundary_xy(idomain_clipped > 0) + mask_topsystem(clipped, ~active_bounds_clipped) + clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty) return clipped diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py index 9cd2663fa..54482791b 100644 --- a/imod/tests/test_mf6/test_ex01_twri.py +++ b/imod/tests/test_mf6/test_ex01_twri.py @@ -577,6 +577,17 @@ def test_slice_and_run_with_state(transient_twri_model_extended, tmp_path): np_array = clipped_boundary["head"].values assert (np_array == 1.23).sum() == 33 + # Test that topsystem packages are masked. + topsystem_keys = ["rch", "drn"] + topsystem_mask = clipped_boundary["head"].notnull().compute() + for key in topsystem_keys: + topsystem_pkg = clipped_simulation["GWF_1"][key] + # Take first grid var + gridded_var = topsystem_pkg.dataset[topsystem_pkg._period_data[0]].compute() + # True wherever topsystem is masked but gridded_var still has a value + bad = gridded_var.notnull() & topsystem_mask + assert not bad.any().item() + @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7") def test_slice_and_run_purge_empty_package(transient_twri_model, tmp_path): From 433cd1000bc93a64dc72fcd14d983378678996ef Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 14:10:42 +0200 Subject: [PATCH 09/39] Add test for masking the topsystem --- .../test_mf6/test_utilities/test_mask.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 imod/tests/test_mf6/test_utilities/test_mask.py diff --git a/imod/tests/test_mf6/test_utilities/test_mask.py b/imod/tests/test_mf6/test_utilities/test_mask.py new file mode 100644 index 000000000..19275469a --- /dev/null +++ b/imod/tests/test_mf6/test_utilities/test_mask.py @@ -0,0 +1,19 @@ +from imod.mf6.utilities.mask import mask_topsystem +from imod.typing.grid import zeros_like + + +def test_mask_topsystem(twri_model): + """ + Test the mask_topsystem utility function by deactivating all cells in the + grid. + """ + # Arrange + gwf_model = twri_model["GWF_1"] + mask = zeros_like(gwf_model.domain) + # Act + mask_topsystem(gwf_model, mask) + # Assert + for key in ["rch", "drn"]: + pkg = gwf_model[key] + gridded_var = pkg.dataset[pkg._period_data[0]].compute() + assert not gridded_var.notnull().any().item() From 32e6772ae0023214e5733f6faf04080d61fca09d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 14:13:13 +0200 Subject: [PATCH 10/39] Return None in type annotation --- imod/mf6/utilities/mask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py index 574422d48..3de4f422d 100644 --- a/imod/mf6/utilities/mask.py +++ b/imod/mf6/utilities/mask.py @@ -3,7 +3,7 @@ from imod.typing import GridDataArray -def mask_topsystem(model: IModel, is_active: GridDataArray): +def mask_topsystem(model: IModel, is_active: GridDataArray) -> None: """ Mask all top system packages in the model inplace with a boolean mask indicating active cells. From b3a59fc32140863a541ab27e4f0361a7f8056bf6 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 14:25:14 +0200 Subject: [PATCH 11/39] Use mask_packages method --- imod/mf6/utilities/mask.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py index 3de4f422d..99a969d9d 100644 --- a/imod/mf6/utilities/mask.py +++ b/imod/mf6/utilities/mask.py @@ -21,5 +21,4 @@ def mask_topsystem(model: IModel, is_active: GridDataArray) -> None: for key, pkg in model.items() if isinstance(pkg, ITopSystemBoundaryCondition) ] - for key in topsystem_packages: - model[key] = model[key].mask(is_active) + model.mask_packages(topsystem_packages, is_active) From 698f36e43ca45fc466dfb7a23420b26d766eab01 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 15:07:17 +0200 Subject: [PATCH 12/39] Regrid iMOD5 IBOUND data when regridding cap data and also mask where ibound is not active --- imod/common/utilities/regrid.py | 4 ++++ imod/mf6/rch.py | 8 +++++--- imod/msw/grid_data.py | 3 ++- imod/msw/regrid/regrid_schemes.py | 1 + imod/msw/utilities/imod5_converter.py | 4 +++- imod/typing/__init__.py | 1 + 6 files changed, 16 insertions(+), 5 deletions(-) diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py index aaf28abb9..f7c389ebf 100644 --- a/imod/common/utilities/regrid.py +++ b/imod/common/utilities/regrid.py @@ -492,9 +492,13 @@ def regrid_imod5_cap_data( cap_data_regridded = _regrid_package_data( imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache ) + bnd_data_regridded = _regrid_package_data( + imod5_data["bnd"], target_grid, regridder_types, regrid_cache + ) extra_paths = imod5_data["extra"]["paths"] imod5_regridded: Imod5DataDict = { "cap": cap_data_regridded, + "bnd": bnd_data_regridded, "extra": {"paths": extra_paths}, } return imod5_regridded diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index 1933f480e..b221b451a 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -298,12 +298,14 @@ def from_imod5_cap_data( used to couple MODFLOW6 to MetaSWAP models. Active cells will have a recharge rate of 0.0. """ - cap_data = regrid_imod5_cap_data( + imod5_data_regridded = regrid_imod5_cap_data( imod5_data, target_dis, regridder_types, regrid_cache - )["cap"] + ) + cap_data = imod5_data_regridded["cap"] + bnd_data = imod5_data_regridded["bnd"] msw_area = get_cell_area_from_imod5_data(cap_data) - msw_active = is_msw_active_cell(target_dis, cap_data, msw_area) + msw_active = is_msw_active_cell(target_dis, cap_data, msw_area, bnd_data) active = msw_active.all data = {} diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index 24fbece43..a15fe8a47 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -202,6 +202,7 @@ def from_imod5_data( as aggregated over subunits. """ imod5_cap = imod5_data["cap"] + imod5_bnd = imod5_data["bnd"] data = {} data["area"] = get_cell_area_from_imod5_data(imod5_cap) @@ -210,7 +211,7 @@ def from_imod5_data( data["surface_elevation"] = imod5_cap["surface_elevation"] data["soil_physical_unit"] = imod5_cap["soil_physical_unit"].astype(int) - msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"], imod5_bnd) data_active = mask_and_broadcast_pkg_data(cls, data, msw_active) data_active["active"] = msw_active.all return cls(**data_active), msw_active diff --git a/imod/msw/regrid/regrid_schemes.py b/imod/msw/regrid/regrid_schemes.py index 30547abfc..a7e73207a 100644 --- a/imod/msw/regrid/regrid_schemes.py +++ b/imod/msw/regrid/regrid_schemes.py @@ -73,6 +73,7 @@ class CapDataRegridMethod(DataclassType): steering_location: RegridVarType = (RegridderType.OVERLAP, "mode") plot_drainage_level: RegridVarType = (RegridderType.OVERLAP, "mean") plot_drainage_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + ibound: RegridVarType = (RegridderType.OVERLAP, "mode") @dataclass(config=_CONFIG) diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 17e2923e3..772f11277 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -96,6 +96,7 @@ def is_msw_active_cell( target_dis: StructuredDiscretization, imod5_cap: GridDataDict, msw_area: GridDataArray, + imod5_bnd: GridDataArray, ) -> MetaSwapActive: """ Return grid of cells that are active in the coupled computation, based on @@ -113,7 +114,8 @@ 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"] > 0) & (msw_area > 0) & (mf6_top_active > 0) + imod5_active = (imod5_bnd["ibound"] > 0) & (imod5_cap["boundary"] > 0) + subunit_active = imod5_active & (msw_area > 0) & (mf6_top_active > 0) active = subunit_active.any(dim="subunit") return MetaSwapActive(active, subunit_active) diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 1c2a4885c..250d9b69a 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -33,6 +33,7 @@ class DropVarsType(TypedDict, total=False): class Imod5DataDict(TypedDict, total=False): + bnd: GridDataDict cap: GridDataDict extra: dict[str, list[str]] From 03aabc0578cd197e1a0f731049632488104e0c37 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 15:09:01 +0200 Subject: [PATCH 13/39] Also drop bnd layer --- imod/util/dims.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/imod/util/dims.py b/imod/util/dims.py index c9228ee39..de3c1fd01 100644 --- a/imod/util/dims.py +++ b/imod/util/dims.py @@ -57,4 +57,8 @@ def _drop_layer_if_dataarray( def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict: cap_data = imod5_data["cap"] - return {"cap": {key: _drop_layer_if_dataarray(da) for key, da in cap_data.items()}} + bnd_data = imod5_data["bnd"] + return { + "cap": {key: _drop_layer_if_dataarray(da) for key, da in cap_data.items()}, + "bnd": {key: _drop_layer_if_dataarray(da) for key, da in bnd_data.items()}, + } From 3171d54129c34848b20ad54b848ce5bd679cd409 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 16:02:37 +0200 Subject: [PATCH 14/39] Call correct var an improve varname --- imod/common/utilities/regrid.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py index f7c389ebf..fd264378a 100644 --- a/imod/common/utilities/regrid.py +++ b/imod/common/utilities/regrid.py @@ -486,14 +486,14 @@ def regrid_imod5_cap_data( and ``imod.mf6.Recharge.from_imod5_cap_data``. """ # Drop layer coords - imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data) + imod5_no_layer = drop_layer_dim_cap_data(imod5_data) target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True) # Regrid the input data cap_data_regridded = _regrid_package_data( - imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache + imod5_no_layer["cap"], target_grid, regridder_types, regrid_cache ) bnd_data_regridded = _regrid_package_data( - imod5_data["bnd"], target_grid, regridder_types, regrid_cache + imod5_no_layer["bnd"], target_grid, regridder_types, regrid_cache ) extra_paths = imod5_data["extra"]["paths"] imod5_regridded: Imod5DataDict = { From ac10b19cce1f7ec478f39a0febed8b0d44d8cd75 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 16:02:46 +0200 Subject: [PATCH 15/39] Add docstring --- imod/msw/utilities/imod5_converter.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 772f11277..7f18a4acf 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -114,6 +114,8 @@ def is_msw_active_cell( Cells active per subunit """ mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True) + # Where IBOUND = -1, there also shouldn't be any active cells in the CAP + # boundary array. imod5_active = (imod5_bnd["ibound"] > 0) & (imod5_cap["boundary"] > 0) subunit_active = imod5_active & (msw_area > 0) & (mf6_top_active > 0) active = subunit_active.any(dim="subunit") From 65cd181622298d2a859a5f9844f7db1dee42b061 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 16:03:41 +0200 Subject: [PATCH 16/39] Add ibound to test fixture and expand tests to test for cell inactivity --- imod/tests/fixtures/msw_imod5_cap_fixture.py | 13 +++++++++++++ imod/tests/test_msw/test_grid_data.py | 17 +++++++++++++++-- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py index 28c53335b..76bb895f6 100644 --- a/imod/tests/fixtures/msw_imod5_cap_fixture.py +++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py @@ -232,6 +232,19 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) + d2 = {} + d2["ibound"] = xr.DataArray( + np.array([ + [ + [1, 1, 1], + [0, 0, 0], + [1, 1, 0], + ]], + dtype=int), + **da_kwargs + ) # fmt: on imod5_data["cap"] = d + imod5_data["bnd"] = d2 + return imod5_data diff --git a/imod/tests/test_msw/test_grid_data.py b/imod/tests/test_msw/test_grid_data.py index f38559d14..c3b40d657 100644 --- a/imod/tests/test_msw/test_grid_data.py +++ b/imod/tests/test_msw/test_grid_data.py @@ -474,7 +474,12 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]): cap_data["soil_physical_unit"] = xr.ones_like(like, dtype=int) cap_data["active"] = xr.ones_like(like, dtype=bool) - imod5_data = {"cap": cap_data} + ibound = xr.ones_like(like, dtype=int) + ibound[0, 0] = -1 + + bnd_data = {} + bnd_data["ibound"] = ibound + imod5_data = {"cap": cap_data, "bnd": bnd_data} layer = xr.DataArray([1, 1], coords={"layer": [1, 2]}, dims=("layer",)) idomain = layer * xr.ones_like(like, dtype=int) @@ -484,7 +489,15 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]): griddata, _ = GridData.from_imod5_data(imod5_data, target_dis=dis) expected_rootzone_depth = cap_data["rootzone_thickness"] * 0.01 + expected_rootzone_depth[0, 0] = np.nan xr.testing.assert_allclose( expected_rootzone_depth, griddata["rootzone_depth"].sel(subunit=0, drop=True) ) - assert (griddata["landuse"].sel(subunit=1, drop=True) == 18).all() + # Test if all cells in subunit = 1 set to "urban" landuse code + np.testing.assert_array_equal( + np.unique(griddata["landuse"].sel(subunit=1, drop=True)), [0, 18] + ) + # Test if cell where IBOUND == -1 is set to inactive (landuse = 0) + np.testing.assert_array_equal(griddata["landuse"][:, 0, 0], [0, 0]) + np.testing.assert_array_equal(griddata["rootzone_depth"][:, 0, 0], [np.nan, np.nan]) + np.testing.assert_array_equal(griddata["active"][0, 0], [False, False]) From 9441fd40c915eaf25520d1ab15b648887a45f8c3 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 16:14:24 +0200 Subject: [PATCH 17/39] Update changelog --- docs/api/changelog.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 6f02e203c..1ecc7c0d5 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -29,6 +29,12 @@ Fixed - Fixed resampling in :meth:`imod.mf6.Well.from_imod5_data` and :meth:`imod.mf6.LayeredWell.from_imod5_data` when simulation timesteps precede the first well timestep. +- :meth:`imod.mf6.GroundwaterFlowModel.from_imod5_data` now masks cells in + topsystem packages (:class:`imod.mf6.River`, + :class:`imod.mf6.GeneralHeadBoundary`, :class:`imod.mf6.Drainage`, + :class:`imod.mf6.Recharge`) where IBOUND is less than 0. +- :meth:`imod.msw.MetaSwapModel.from_imod5_data` now masks cells where IBOUND is + less than 0. Changed ~~~~~~~ @@ -36,6 +42,11 @@ Changed - Deprecated :class:`imod.msw.Sprinkling` in favor of :class:`imod.msw.SprinklingGrid`. Call :class:`imod.msw.SprinklingGrid` to get the same behavior as you were used to. +- If ``states_for_boundary`` is provided to + :meth:`imod.mf6.GroundwaterFlowModel.clip_box`, topsystem packages + (:class:`imod.mf6.River`, :class:`imod.mf6.GeneralHeadBoundary`, + :class:`imod.mf6.Drainage`, :class:`imod.mf6.Recharge`) will also be masked + where constant head cells are placed. [1.1.0] - 2026-08-03 -------------------- From e6e7a0365a21d1b7c09ff538ba0ea22244e0f08d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 17:00:00 +0200 Subject: [PATCH 18/39] Rename to avoid duplicate test module names --- .../test_utilities/{test_mask.py => test_mf6_mask_util.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename imod/tests/test_mf6/test_utilities/{test_mask.py => test_mf6_mask_util.py} (100%) diff --git a/imod/tests/test_mf6/test_utilities/test_mask.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py similarity index 100% rename from imod/tests/test_mf6/test_utilities/test_mask.py rename to imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py From fb1b855660fb95d5f2e98dc13d74c298d955721b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 17:03:16 +0200 Subject: [PATCH 19/39] Also rename msw mask util test module --- .../test_utilities/{test_mask.py => test_msw_mask_util.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename imod/tests/test_msw/test_utilities/{test_mask.py => test_msw_mask_util.py} (100%) diff --git a/imod/tests/test_msw/test_utilities/test_mask.py b/imod/tests/test_msw/test_utilities/test_msw_mask_util.py similarity index 100% rename from imod/tests/test_msw/test_utilities/test_mask.py rename to imod/tests/test_msw/test_utilities/test_msw_mask_util.py From 52b6ea42c5af451ce3c8b9000dd340559d616c62 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 17:18:33 +0200 Subject: [PATCH 20/39] Fix and expand mf6 mask tests --- .../test_utilities/test_mf6_mask_util.py | 31 ++++++++++++++++--- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py index 19275469a..8993fa7b3 100644 --- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py +++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py @@ -1,6 +1,6 @@ from imod.mf6.utilities.mask import mask_topsystem -from imod.typing.grid import zeros_like - +from imod.typing.grid import zeros_like, ones_like +import numpy as np def test_mask_topsystem(twri_model): """ @@ -9,11 +9,32 @@ def test_mask_topsystem(twri_model): """ # Arrange gwf_model = twri_model["GWF_1"] - mask = zeros_like(gwf_model.domain) + is_active = ones_like(gwf_model.domain) + # Mask first cell + is_active[0, 0, 0] = 0 + # Act - mask_topsystem(gwf_model, mask) + mask_topsystem(gwf_model, is_active) # Assert for key in ["rch", "drn"]: pkg = gwf_model[key] gridded_var = pkg.dataset[pkg._period_data[0]].compute() - assert not gridded_var.notnull().any().item() + first_cell = gridded_var.data.ravel()[0] + assert np.isnan(first_cell).item() + + + +def test_mask_topsystem__all_removed(twri_model): + """ + Test the mask_topsystem utility function by deactivating all cells in the + grid. + """ + # Arrange + gwf_model = twri_model["GWF_1"] + is_active = zeros_like(gwf_model.domain) + # Act + mask_topsystem(gwf_model, is_active) + # Assert + for key in ["rch", "drn"]: + assert key not in gwf_model.keys() + From fc7d39dcc4659a86f6b7a1bd210034355046e35b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 17:27:35 +0200 Subject: [PATCH 21/39] Update mock setup --- imod/tests/test_mf6/test_mf6_model.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py index b7a3ecd4c..cea2b47ee 100644 --- a/imod/tests/test_mf6/test_mf6_model.py +++ b/imod/tests/test_mf6/test_mf6_model.py @@ -263,9 +263,20 @@ def test_clip_box_with_state_for_boundary( # Arrange. state_for_boundary = MagicMock(spec_set=UgridDataArray) + idomain = xr.DataArray( + np.ones((1, 2, 2), dtype=np.int32), dims=("layer", "y", "x") + ) + top = xr.DataArray(np.ones((2, 2), dtype=np.float64), dims=("y", "x")) + bottom = xr.DataArray(np.array([-1.0], dtype=np.float64), dims=("layer",)) + discretization_mock = MagicMock(spec_set=Package) discretization_mock._pkg_id = "dis" discretization_mock.clip_box.return_value = discretization_mock + discretization_mock.__getitem__.side_effect = { + "idomain": idomain, + "top": top, + "bottom": bottom, + }.__getitem__ clipped_boundary_mock = MagicMock(spec_set=pkg_type) clipped_boundary_mock.is_empty.return_value = False @@ -305,10 +316,21 @@ def test_clip_box_with_unassigned_boundaries_in_original_model( # Arrange. state_for_boundary = MagicMock(spec_set=UgridDataArray) + idomain = xr.DataArray( + np.ones((1, 2, 2), dtype=np.int32), dims=("layer", "y", "x") + ) + top = xr.DataArray(np.ones((2, 2), dtype=np.float64), dims=("y", "x")) + bottom = xr.DataArray(np.array([-1.0], dtype=np.float64), dims=("layer",)) + discretization_mock = MagicMock(spec_set=Package) discretization_mock._pkg_id = "dis" discretization_mock.is_empty.side_effect = [False, False] discretization_mock.clip_box.return_value = discretization_mock + discretization_mock.__getitem__.side_effect = { + "idomain": idomain, + "top": top, + "bottom": bottom, + }.__getitem__ constant_boundary_mock = MagicMock(spec_set=pkg_type) constant_boundary_mock.is_empty.side_effect = [False, False] From 4e060446d94f8390b43b27f0502ec64fac319712 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 23 Sep 2026 17:28:13 +0200 Subject: [PATCH 22/39] format --- imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py index 8993fa7b3..eded019f6 100644 --- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py +++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py @@ -1,7 +1,9 @@ -from imod.mf6.utilities.mask import mask_topsystem -from imod.typing.grid import zeros_like, ones_like import numpy as np +from imod.mf6.utilities.mask import mask_topsystem +from imod.typing.grid import ones_like, zeros_like + + def test_mask_topsystem(twri_model): """ Test the mask_topsystem utility function by deactivating all cells in the @@ -23,7 +25,6 @@ def test_mask_topsystem(twri_model): assert np.isnan(first_cell).item() - def test_mask_topsystem__all_removed(twri_model): """ Test the mask_topsystem utility function by deactivating all cells in the @@ -37,4 +38,3 @@ def test_mask_topsystem__all_removed(twri_model): # Assert for key in ["rch", "drn"]: assert key not in gwf_model.keys() - From cb51c400c64bc6dec428d39a671b2d1e02dc8229 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 11:51:58 +0200 Subject: [PATCH 23/39] Add ibound to regrid schemes where it was missing and slightly improve docstring --- imod/mf6/regrid/regrid_schemes.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/imod/mf6/regrid/regrid_schemes.py b/imod/mf6/regrid/regrid_schemes.py index f2a957316..1348488ff 100644 --- a/imod/mf6/regrid/regrid_schemes.py +++ b/imod/mf6/regrid/regrid_schemes.py @@ -456,13 +456,15 @@ class StorageCoefficientRegridMethod(DataclassType): class CapDataRechargeRegridMethod(DataclassType): """ Object containing regridder methods for CAP data for the - :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains regridder - methods for only the relevant CAP variables for the recharge package. + :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains + regridder methods for only the relevant iMOD5 CAP and BND variables for the + recharge package. """ boundary: RegridVarType = (RegridderType.OVERLAP, "mode") wetted_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") urban_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") + ibound: RegridVarType = (RegridderType.OVERLAP, "mode") @dataclass(config=_CONFIG) @@ -470,8 +472,10 @@ class CapDataWellRegridMethod(DataclassType): """ Object containing regridder methods for CAP data for the :class:`imod.mf6.LayeredWell.from_imod5_cap_data` method. This contains - regridder methods for only the relevant CAP variables for the well package. + regridder methods for only the relevant iMOD5 CAP and BND variables for the + well package. """ artificial_recharge: RegridVarType = (RegridderType.OVERLAP, "mean") artificial_recharge_layer: RegridVarType = (RegridderType.OVERLAP, "mode") + ibound: RegridVarType = (RegridderType.OVERLAP, "mode") From 4776d41dfde075654c7327f09576f7f09b0c8cb4 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 11:52:18 +0200 Subject: [PATCH 24/39] Include bnd ibound data in test fixtures where missing. --- imod/tests/fixtures/imod5_cap_data.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py index 8e9b29a9b..3330a4099 100644 --- a/imod/tests/fixtures/imod5_cap_data.py +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -82,8 +82,8 @@ def cap_data_sprinkling_grid() -> Imod5DataDict: "artificial_recharge_layer": layer, "artificial_recharge_capacity": xr.DataArray(25.0), } - - return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} + bnd_data = {"ibound": zeros_grid(n) + 1} + return {"cap": cap_data, "bnd": bnd_data, "extra": {"paths": ["path1", "path2"]}} @pytest.fixture(scope="function") @@ -105,8 +105,9 @@ def cap_data_sprinkling_grid__big() -> Imod5DataDict: "artificial_recharge_layer": layer, "artificial_recharge_capacity": xr.DataArray(25.0), } + bnd_data = {"ibound": zeros_dask_grid(n) + 1} - return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} + return {"cap": cap_data, "bnd": bnd_data, "extra": {"paths": ["path1", "path2"]}} @pytest.fixture(scope="function") From 65ec47d31a2acfc76fc3a70c8f4f06aa5b3b9c0a Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 14:52:52 +0200 Subject: [PATCH 25/39] Include ignore_time_purge_empty in mask_topsystem calls. Fix creation of mask for topsystem in clip_box. --- imod/mf6/model.py | 25 +++++++++++++++---------- imod/mf6/model_gwf.py | 1 + imod/mf6/utilities/imod5_converter.py | 3 ++- imod/mf6/utilities/mask.py | 8 ++++++-- 4 files changed, 24 insertions(+), 13 deletions(-) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 70f21325f..4f3400665 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -49,7 +49,6 @@ from imod.mf6.wel import GridAgnosticWell from imod.mf6.write_context import WriteContext from imod.schemata import SchemataDict, ValidationError -from imod.select.grid import active_grid_boundary_xy from imod.typing import GridDataArray from imod.util.regrid import RegridderWeightsCache @@ -810,18 +809,24 @@ def clip_box( clipped_boundary_condition = _create_boundary_condition_clipped_boundary( self, clipped, state_for_boundary, clip_box_args ) - state_pkg_id = self._boundary_state_pkg_type._pkg_id - pkg_name = f"{state_pkg_id}_clipped" if clipped_boundary_condition is not None: - clipped[pkg_name] = clipped_boundary_condition + # Assign clipped boundary condition package + state_pkg_id = self._boundary_state_pkg_type._pkg_id + pkg_name = f"{state_pkg_id}_clipped" - # Clip topsystem packages where the active grid boundary has - # changed. - _, _, idomain_clipped = clipped._get_domain_geometry() - active_bounds_clipped = active_grid_boundary_xy(idomain_clipped > 0) - mask_topsystem(clipped, ~active_bounds_clipped) + clipped[pkg_name] = clipped_boundary_condition - clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty) + # Mask topsystem packages where the state boundary cells have been + # added. + state_varname = clipped_boundary_condition._period_data[0] + state_var = clipped_boundary_condition.dataset[state_varname].isel( + time=0, missing_dims="ignore" + ) + not_added_bc = np.isnan(state_var) + # Purge empty packages called by the mask_topsystem function + mask_topsystem(clipped, not_added_bc, ignore_time_purge_empty) + else: + clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty) return clipped diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 5296bdc62..906d1e4bb 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -460,5 +460,6 @@ def from_imod5_data( result, cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")), regrid_cache, + ignore_time_purge_empty=True, ) return result diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index 0da531a2a..81c5e84d7 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -192,6 +192,7 @@ def mask_topsystem_packages_with_ibound( model: IModel, regridder_types: Optional[ConstantHeadRegridMethod], regrid_cache: RegridderWeightsCache, + ignore_time_purge_empty: bool, ) -> None: """ Mask all top system packages where IBOUND < 0. These locations are assigned @@ -211,4 +212,4 @@ def mask_topsystem_packages_with_ibound( )["ibound"] is_active = regridded_ibound >= 0 - mask_topsystem(model, is_active) + mask_topsystem(model, is_active, ignore_time_purge_empty) diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py index 99a969d9d..d734d5518 100644 --- a/imod/mf6/utilities/mask.py +++ b/imod/mf6/utilities/mask.py @@ -3,7 +3,9 @@ from imod.typing import GridDataArray -def mask_topsystem(model: IModel, is_active: GridDataArray) -> None: +def mask_topsystem( + model: IModel, is_active: GridDataArray, ignore_time_purge_empty: bool +) -> None: """ Mask all top system packages in the model inplace with a boolean mask indicating active cells. @@ -15,10 +17,12 @@ def mask_topsystem(model: IModel, is_active: GridDataArray) -> None: is_active : GridDataArray A boolean array indicating active cells. Top system packages will be masked where this array is False. + ignore_time_purge_empty : bool + If True, ignore the time dimension when masking the packages. """ topsystem_packages = [ key for key, pkg in model.items() if isinstance(pkg, ITopSystemBoundaryCondition) ] - model.mask_packages(topsystem_packages, is_active) + model.mask_packages(topsystem_packages, is_active, ignore_time_purge_empty) From 50f71e5f5965e7a8950b99ccc9f245ea80a435f8 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 15:16:39 +0200 Subject: [PATCH 26/39] Refactor: Move boundary condition creator utiltity functions from model module to their respective --- imod/common/interfaces/imodel.py | 6 + imod/mf6/model.py | 90 +------------- imod/mf6/utilities/clipped_bc_creator.py | 113 +++++++++++++++++- ..._mf6_clipped_boundary_condition_creator.py | 6 +- 4 files changed, 122 insertions(+), 93 deletions(-) diff --git a/imod/common/interfaces/imodel.py b/imod/common/interfaces/imodel.py index 0b30e0a6e..fd143feb7 100644 --- a/imod/common/interfaces/imodel.py +++ b/imod/common/interfaces/imodel.py @@ -2,6 +2,7 @@ from typing import Any, Optional, Tuple from imod.common.interfaces.idict import IDict +from imod.common.interfaces.ipackage import IPackage from imod.common.statusinfo import StatusInfoBase from imod.mf6.validation_settings import ValidationSettings from imod.typing import GridDataArray @@ -68,3 +69,8 @@ def _is_splitting_supported(self) -> Tuple[bool, str]: @abstractmethod def _is_clipping_supported(self) -> Tuple[bool, str]: raise NotImplementedError + + @property + @abstractmethod + def _boundary_state_pkg_type(self) -> type[IPackage]: + raise NotImplementedError diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 4f3400665..83bfa725d 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -20,7 +20,6 @@ from imod.common.interfaces.imodel import IModel from imod.common.serializer import EngineType from imod.common.statusinfo import NestedStatusInfo, StatusInfo, StatusInfoBase -from imod.common.utilities.clip import clip_box_dataset from imod.common.utilities.dump_model import dump_model from imod.common.utilities.mask import mask_packages from imod.common.utilities.regrid import _regrid_like @@ -40,8 +39,7 @@ from imod.mf6.riv import River from imod.mf6.utilities.clipped_bc_creator import ( StateClassType, - StateType, - create_clipped_boundary, + create_boundary_condition_clipped_boundary, ) from imod.mf6.utilities.mask import mask_topsystem from imod.mf6.utilities.mf6hfb import merge_hfb_packages @@ -63,90 +61,6 @@ def pkg_has_cleanup(pkg: Package): return any(isinstance(pkg, pkgtype) for pkgtype in PKGTYPES_WITH_CLEANUP) -def _create_boundary_condition_for_unassigned_boundary( - model: Modflow6Model, - state_for_boundary: Optional[GridDataArray], - additional_boundaries: list[Optional[StateType]] = [None], -) -> Optional[StateType]: - if state_for_boundary is None: - return None - - pkg_type = model._boundary_state_pkg_type - constant_state_packages = [ - pkg for _, pkg in model.items() if isinstance(pkg, pkg_type) - ] - - filtered_boundaries: list[StateType] = [ - item for item in additional_boundaries or [] if item is not None - ] - - constant_state_packages.extend(filtered_boundaries) - - return create_clipped_boundary( - model.domain, state_for_boundary, constant_state_packages, pkg_type - ) - - -def _create_boundary_condition_clipped_boundary( - original_model: Modflow6Model, - clipped_model: Modflow6Model, - state_for_boundary: Optional[GridDataArray], - clip_box_args: tuple[Any, ...], -) -> Optional[StateType]: - # Create temporary boundary condition for the original model boundary. This - # is used later to see which boundaries can be ignored as they were already - # present in the original model. We want to just end up with the boundary - # created by the clip. - unassigned_boundary_original_domain = ( - _create_boundary_condition_for_unassigned_boundary( - original_model, state_for_boundary - ) - ) - # Clip the unassigned boundary to the clipped model's domain, required to - # avoid topological errors later. - if unassigned_boundary_original_domain is not None: - unassigned_boundary_clipped = unassigned_boundary_original_domain.clip_box( - *clip_box_args - ) - else: - unassigned_boundary_clipped = None - - if state_for_boundary is not None: - # Clip box as dataset, temporarily add variable name to convert to - # dataset, then turn back into DataArray. - varname = original_model._boundary_state_pkg_type._period_data[0] - state_for_boundary = state_for_boundary.to_dataset(name=varname) - state_for_boundary_clipped = clip_box_dataset( - state_for_boundary, *clip_box_args - )[varname] - else: - state_for_boundary_clipped = None - - bc_constant_pkg = _create_boundary_condition_for_unassigned_boundary( - clipped_model, state_for_boundary_clipped, [unassigned_boundary_clipped] - ) - - # Remove all indices before first timestep of state_for_clipped_boundary. - # This to prevent empty dataarrays unnecessarily being made for these - # indices, which can lead to them to be removed when purging empty packages - # with ignore_time=True. Unfortunately, this is needs to be handled here and - # not in _create_boundary_condition_for_unassigned_boundary, as otherwise - # this function is called twice which could result in broadcasting errors in - # the second call if the time domain of state_for_boundary and assigned - # packages have no overlap. - if ( - (state_for_boundary is not None) - and (state_for_boundary.indexes.get("time") is not None) - and (bc_constant_pkg is not None) - ): - start_time = state_for_boundary.indexes["time"][0] - bc_constant_pkg.dataset = bc_constant_pkg.dataset.sel( - time=slice(start_time, None) - ) - - return bc_constant_pkg - - class Modflow6Model(collections.UserDict[str, Package], IModel, abc.ABC): _mandatory_packages: tuple[str, ...] = () _init_schemata: SchemataDict = {} @@ -806,7 +720,7 @@ def clip_box( *clip_box_args, ) - clipped_boundary_condition = _create_boundary_condition_clipped_boundary( + clipped_boundary_condition = create_boundary_condition_clipped_boundary( self, clipped, state_for_boundary, clip_box_args ) if clipped_boundary_condition is not None: diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py index ea6f0f415..c666ae261 100644 --- a/imod/mf6/utilities/clipped_bc_creator.py +++ b/imod/mf6/utilities/clipped_bc_creator.py @@ -1,7 +1,9 @@ -from typing import Optional, Tuple, TypeAlias +from typing import Any, Optional, Tuple, TypeAlias, cast import xarray as xr +from imod.common.interfaces.imodel import IModel +from imod.common.utilities.clip import clip_box_dataset from imod.mf6 import ConstantConcentration, ConstantHead from imod.select.grid import active_grid_boundary_xy from imod.typing import GridDataArray @@ -117,7 +119,7 @@ def _create_clipped_boundary_state( return state_for_clipped_boundary.where(unassigned_grid_boundaries) -def create_clipped_boundary( +def _create_clipped_boundary_pkg( idomain: GridDataArray, state_for_clipped_boundary: GridDataArray, original_constant_head_boundaries: list[StateType], @@ -152,3 +154,110 @@ def create_clipped_boundary( ) return pkg_type(constant_state, print_input=True, print_flows=True, save_flows=True) + + +def _create_boundary_condition_for_unassigned_boundary( + model: IModel, + state_for_boundary: Optional[GridDataArray], + additional_boundaries: list[Optional[StateType]] = [None], +) -> Optional[StateType]: + if state_for_boundary is None: + return None + + pkg_type = cast(StateClassType, model._boundary_state_pkg_type) + constant_state_packages = [ + pkg for _, pkg in model.items() if isinstance(pkg, pkg_type) + ] + + filtered_boundaries: list[StateType] = [ + item for item in additional_boundaries or [] if item is not None + ] + + constant_state_packages.extend(filtered_boundaries) + + return _create_clipped_boundary_pkg( + model.domain, state_for_boundary, constant_state_packages, pkg_type + ) + + +def create_boundary_condition_clipped_boundary( + original_model: IModel, + clipped_model: IModel, + state_for_boundary: Optional[GridDataArray], + clip_box_args: tuple[Any, ...], +) -> Optional[StateType]: + """ + Create a clipped boundary condition for a given state in the clipped model. + The function takes the original model as a reference to determine where + boundary conditions should NOT be placed, then applies this information to + create the boundary condition in the clipped model. + + Parameters + ---------- + original_model : IModel + The original model containing the unassigned boundary condition. + clipped_model : IModel + The clipped model where the boundary condition will be applied. + state_for_boundary : Optional[GridDataArray] + The state array for the boundary condition. + clip_box_args : tuple[Any, ...] + Arguments defining the clipping box. + + Returns + ------- + Optional[StateType] + The clipped boundary condition package, or None if no boundary condition is created. + """ + # Create temporary boundary condition for the original model boundary. This + # is used later to see which boundaries can be ignored as they were already + # present in the original model. We want to just end up with the boundary + # created by the clip. + unassigned_boundary_original_domain = ( + _create_boundary_condition_for_unassigned_boundary( + original_model, state_for_boundary + ) + ) + # Clip the unassigned boundary to the clipped model's domain, required to + # avoid topological errors later. + if unassigned_boundary_original_domain is not None: + unassigned_boundary_clipped = unassigned_boundary_original_domain.clip_box( + *clip_box_args + ) + else: + unassigned_boundary_clipped = None + + if state_for_boundary is not None: + # Clip box as dataset, temporarily add variable name to convert to + # dataset, then turn back into DataArray. + state_cls = cast(StateClassType, original_model._boundary_state_pkg_type) + varname = state_cls._period_data[0] + state_for_boundary = state_for_boundary.to_dataset(name=varname) + state_for_boundary_clipped = clip_box_dataset( + state_for_boundary, *clip_box_args + )[varname] + else: + state_for_boundary_clipped = None + + bc_constant_pkg = _create_boundary_condition_for_unassigned_boundary( + clipped_model, state_for_boundary_clipped, [unassigned_boundary_clipped] + ) + + # Remove all indices before first timestep of state_for_clipped_boundary. + # This to prevent empty dataarrays unnecessarily being made for these + # indices, which can lead to them to be removed when purging empty packages + # with ignore_time=True. Unfortunately, this is needs to be handled here and + # not in _create_boundary_condition_for_unassigned_boundary, as otherwise + # this function is called twice which could result in broadcasting errors in + # the second call if the time domain of state_for_boundary and assigned + # packages have no overlap. + if ( + (state_for_boundary is not None) + and (state_for_boundary.indexes.get("time") is not None) + and (bc_constant_pkg is not None) + ): + start_time = state_for_boundary.indexes["time"][0] + bc_constant_pkg.dataset = bc_constant_pkg.dataset.sel( + time=slice(start_time, None) + ) + + return bc_constant_pkg diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py b/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py index aa599ebca..2be514194 100644 --- a/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py +++ b/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py @@ -5,7 +5,7 @@ from imod.mf6 import ConstantHead from imod.mf6.utilities.clipped_bc_creator import ( - create_clipped_boundary, + _create_clipped_boundary_pkg, ) from imod.select.grid import grid_boundary_xy @@ -55,7 +55,7 @@ def test_create_different_n_clipped_cells(self, circle_dis, n_clipped_cells): ) # Act. - constant_head_pkg_clipped_domain = create_clipped_boundary( + constant_head_pkg_clipped_domain = _create_clipped_boundary_pkg( idomain, clipped_boundary_values, [reduced_boundary_constant_head_pkg], @@ -104,7 +104,7 @@ def test_create_different_dis(self, dis, grid_data_array, request): ) # Act. - constant_head_pkg_clipped_domain = create_clipped_boundary( + constant_head_pkg_clipped_domain = _create_clipped_boundary_pkg( idomain, clipped_boundary_values, [reduced_boundary_constant_head_pkg], From d8b191ae8912ed1871afbdc28d6c51627bb0b8fb Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 16:06:09 +0200 Subject: [PATCH 27/39] Update mocking framework --- imod/tests/test_mf6/test_mf6_model.py | 34 ++++++++++++++++++++------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py index cea2b47ee..42f19e85b 100644 --- a/imod/tests/test_mf6/test_mf6_model.py +++ b/imod/tests/test_mf6/test_mf6_model.py @@ -102,12 +102,18 @@ def test_circle_roundtrip(circle_model, tmp_path): roundtrip(circle_model["GWF_1"], tmp_path) +class ConcreteModflow6Model(Modflow6Model): + """Concrete implementation of the abstract Modflow6Model for testing purposes.""" + + _boundary_state_pkg_type = ConstantHead + + class TestModel: def test_write_valid_model_without_error(self, tmpdir_factory): # Arrange. tmp_path = tmpdir_factory.mktemp("TestSimulation") model_name = "Test model" - model = Modflow6Model() + model = ConcreteModflow6Model() # create write context validation_context = ValidationSettings() write_context = WriteContext(tmp_path) @@ -135,7 +141,7 @@ def test_write_without_dis_pkg_return_error(self, tmpdir_factory): # Arrange. tmp_path = tmpdir_factory.mktemp("TestSimulation") model_name = "Test model" - model = Modflow6Model() + model = ConcreteModflow6Model() # create write context validation_context = ValidationSettings() write_context = WriteContext(tmp_path) @@ -158,7 +164,7 @@ def test_write_with_invalid_pkg_returns_error(self, tmpdir_factory): # Arrange. tmp_path = tmpdir_factory.mktemp("TestSimulation") model_name = "Test model" - model = Modflow6Model() + model = ConcreteModflow6Model() # create write context validation_context = ValidationSettings() write_context = WriteContext(tmp_path) @@ -192,7 +198,7 @@ def test_write_with_two_invalid_pkg_returns_two_errors(self, tmpdir_factory): validation_context = ValidationSettings() write_context = WriteContext(simulation_directory=tmp_path) - model = Modflow6Model() + model = ConcreteModflow6Model() discretization_mock = MagicMock(spec_set=Package) discretization_mock._pkg_id = "dis" @@ -249,7 +255,8 @@ def test_clip_box_without_state_for_boundary(self, model_type, pkg_type): pkg_id = pkg_type._pkg_id assert f"{pkg_id}_clipped" not in clipped - @mock.patch("imod.mf6.model.create_clipped_boundary") + @mock.patch("imod.mf6.model.mask_topsystem") + @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg") @pytest.mark.parametrize( "model_type, pkg_type", [ @@ -258,7 +265,11 @@ def test_clip_box_without_state_for_boundary(self, model_type, pkg_type): ], ) def test_clip_box_with_state_for_boundary( - self, create_clipped_boundary_mock, model_type, pkg_type + self, + create_clipped_boundary_mock, + mask_topsystem_mock, + model_type, + pkg_type, ): # Arrange. state_for_boundary = MagicMock(spec_set=UgridDataArray) @@ -301,8 +312,10 @@ def test_clip_box_with_state_for_boundary( [], pkg_type, ) + mask_topsystem_mock.assert_called_once() - @mock.patch("imod.mf6.model.create_clipped_boundary") + @mock.patch("imod.mf6.model.mask_topsystem") + @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg") @pytest.mark.parametrize( "model_type, pkg_type", [ @@ -311,7 +324,11 @@ def test_clip_box_with_state_for_boundary( ], ) def test_clip_box_with_unassigned_boundaries_in_original_model( - self, create_clipped_boundary_mock, model_type, pkg_type + self, + create_clipped_boundary_mock, + mask_topsystem_mock, + model_type, + pkg_type, ): # Arrange. state_for_boundary = MagicMock(spec_set=UgridDataArray) @@ -361,6 +378,7 @@ def test_clip_box_with_unassigned_boundaries_in_original_model( [constant_boundary_mock, unassigned_original_constant_boundary.clip_box()], pkg_type, ) + mask_topsystem_mock.assert_called_once() class TestGroundwaterFlowModel: From 0f41608e526be423db4eb7e5570132d4d2d8ccd5 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 16:06:21 +0200 Subject: [PATCH 28/39] Update missing args --- imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py index eded019f6..167d6d1d4 100644 --- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py +++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py @@ -16,7 +16,7 @@ def test_mask_topsystem(twri_model): is_active[0, 0, 0] = 0 # Act - mask_topsystem(gwf_model, is_active) + mask_topsystem(gwf_model, is_active, True) # Assert for key in ["rch", "drn"]: pkg = gwf_model[key] @@ -34,7 +34,7 @@ def test_mask_topsystem__all_removed(twri_model): gwf_model = twri_model["GWF_1"] is_active = zeros_like(gwf_model.domain) # Act - mask_topsystem(gwf_model, is_active) + mask_topsystem(gwf_model, is_active, True) # Assert for key in ["rch", "drn"]: assert key not in gwf_model.keys() From 497f94a55b4563df15667732e74ad0e29b4372bd Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 24 Sep 2026 16:12:17 +0200 Subject: [PATCH 29/39] Drop time coord and add docstring --- imod/mf6/model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 83bfa725d..6d783fd40 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -733,8 +733,10 @@ def clip_box( # Mask topsystem packages where the state boundary cells have been # added. state_varname = clipped_boundary_condition._period_data[0] + # Select the state variable for the first time step as mask. + # Its location will be constant through time. state_var = clipped_boundary_condition.dataset[state_varname].isel( - time=0, missing_dims="ignore" + time=0, missing_dims="ignore", drop=True ) not_added_bc = np.isnan(state_var) # Purge empty packages called by the mask_topsystem function From 0bede8995b1170365440a90c833c79d37c5e4d9f Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 13:13:18 +0200 Subject: [PATCH 30/39] Rename to "regrid_imod5_cap_and_bnd_data" --- imod/common/utilities/regrid.py | 2 +- imod/mf6/rch.py | 4 ++-- imod/mf6/utilities/imod5_converter.py | 7 +++++-- imod/msw/model.py | 4 ++-- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py index fd264378a..7fd220c6a 100644 --- a/imod/common/utilities/regrid.py +++ b/imod/common/utilities/regrid.py @@ -471,7 +471,7 @@ def _get_regridding_domain( return new_idomain -def regrid_imod5_cap_data( +def regrid_imod5_cap_and_bnd_data( imod5_data: Imod5DataDict, target_dis: IRegridPackage, regridder_types: DataclassType, diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index b221b451a..731676776 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -6,7 +6,7 @@ from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType -from imod.common.utilities.regrid import regrid_imod5_cap_data +from imod.common.utilities.regrid import regrid_imod5_cap_and_bnd_data from imod.logging import init_log_decorator from imod.mf6.aggregate.aggregate_schemes import RechargeAggregationMethod from imod.mf6.dis import StructuredDiscretization, VerticesDiscretization @@ -298,7 +298,7 @@ def from_imod5_cap_data( used to couple MODFLOW6 to MetaSWAP models. Active cells will have a recharge rate of 0.0. """ - imod5_data_regridded = regrid_imod5_cap_data( + imod5_data_regridded = regrid_imod5_cap_and_bnd_data( imod5_data, target_dis, regridder_types, regrid_cache ) cap_data = imod5_data_regridded["cap"] diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index 81c5e84d7..ef2f1d20f 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -7,7 +7,10 @@ from imod.common.interfaces.imodel import IModel from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType -from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data +from imod.common.utilities.regrid import ( + _regrid_package_data, + regrid_imod5_cap_and_bnd_data, +) from imod.mf6.package import Package from imod.mf6.regrid.regrid_schemes import ConstantHeadRegridMethod from imod.mf6.utilities.mask import mask_topsystem @@ -136,7 +139,7 @@ def well_from_imod5_cap_data( "target_dis must be provided when converting iMOD5 cap data " "from grids (IDF)" ) - cap_data_regridded = regrid_imod5_cap_data( + cap_data_regridded = regrid_imod5_cap_and_bnd_data( imod5_data, target_dis, regridder_types, regrid_cache )["cap"] return _well_from_imod5_cap_grid_data(cap_data_regridded) diff --git a/imod/msw/model.py b/imod/msw/model.py index 2a4081296..19832db6e 100644 --- a/imod/msw/model.py +++ b/imod/msw/model.py @@ -19,7 +19,7 @@ from imod.common.utilities.clip import clip_by_grid from imod.common.utilities.dump_model import dump_model from imod.common.utilities.partitioninfo import create_partition_info -from imod.common.utilities.regrid import regrid_imod5_cap_data +from imod.common.utilities.regrid import regrid_imod5_cap_and_bnd_data from imod.common.utilities.version import prepend_content_with_version_info from imod.mf6.dis import StructuredDiscretization from imod.mf6.mf6_wel_adapter import Mf6Wel @@ -823,7 +823,7 @@ def from_imod5_data( parasim_settings = read_para_sim(path_to_parasim) unsa_svat_path = cast(str, parasim_settings["unsa_svat_path"]) # Regrid iMOD5 CAP data to target discretization. - imod5_regridded = regrid_imod5_cap_data( + imod5_regridded = regrid_imod5_cap_and_bnd_data( imod5_data, target_dis, regridder_types, regrid_cache ) # Test with regridded data instead of masked, as masking broadcasts From 9fbc7b69b07fdf74acb0b44f976247b2bb3e514f Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 13:53:19 +0200 Subject: [PATCH 31/39] Rename to drop_layer_dim_cap_and_bnd_data --- imod/common/utilities/regrid.py | 4 ++-- imod/util/dims.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py index 7fd220c6a..94ab76f75 100644 --- a/imod/common/utilities/regrid.py +++ b/imod/common/utilities/regrid.py @@ -27,7 +27,7 @@ is_unstructured, ones_like, ) -from imod.util.dims import drop_layer_dim_cap_data, enforced_dim_order +from imod.util.dims import drop_layer_dim_cap_and_bnd_data, enforced_dim_order from imod.util.regrid import ( RegridderType, RegridderWeightsCache, @@ -486,7 +486,7 @@ def regrid_imod5_cap_and_bnd_data( and ``imod.mf6.Recharge.from_imod5_cap_data``. """ # Drop layer coords - imod5_no_layer = drop_layer_dim_cap_data(imod5_data) + imod5_no_layer = drop_layer_dim_cap_and_bnd_data(imod5_data) target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True) # Regrid the input data cap_data_regridded = _regrid_package_data( diff --git a/imod/util/dims.py b/imod/util/dims.py index de3c1fd01..535ef6cb1 100644 --- a/imod/util/dims.py +++ b/imod/util/dims.py @@ -55,7 +55,7 @@ def _drop_layer_if_dataarray( return _drop_layer_dim_and_coord(da) -def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict: +def drop_layer_dim_cap_and_bnd_data(imod5_data: Imod5DataDict) -> Imod5DataDict: cap_data = imod5_data["cap"] bnd_data = imod5_data["bnd"] return { From 6aff0f38d83e2be9e1d82483a63c4f235ac8e2d2 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 14:12:06 +0200 Subject: [PATCH 32/39] Add missing docstrings to functions --- imod/mf6/utilities/clipped_bc_creator.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py index c666ae261..744bd0440 100644 --- a/imod/mf6/utilities/clipped_bc_creator.py +++ b/imod/mf6/utilities/clipped_bc_creator.py @@ -17,6 +17,9 @@ def _find_unassigned_grid_boundaries( active_grid_boundary: GridDataArray, boundary_conditions: list[StateType], ) -> GridDataArray: + """ + Find the grid boundaries that have not been assigned any boundary conditions. + """ unassigned_grid_boundaries = active_grid_boundary.copy() for boundary_condition in boundary_conditions: # Fetch variable name from the first boundary condition, can be "head" or @@ -105,7 +108,7 @@ def _create_clipped_boundary_state( idomain: GridDataArray, state_for_clipped_boundary: GridDataArray, original_constant_head_boundaries: list[StateType], -): +) -> GridDataArray: """Helper function to make sure dimension order is enforced""" active_grid_boundary = active_grid_boundary_xy(idomain > 0) unassigned_grid_boundaries = _find_unassigned_grid_boundaries( @@ -161,6 +164,14 @@ def _create_boundary_condition_for_unassigned_boundary( state_for_boundary: Optional[GridDataArray], additional_boundaries: list[Optional[StateType]] = [None], ) -> Optional[StateType]: + """ + Create a boundary condition for the unnassigned boundary cells of a model. + + Constant state packages are collected from the model and any additional + boundaries provided will be added to this list. Cells which are not covered + by any of these packages will have a new boundary condition created for + them. + """ if state_for_boundary is None: return None @@ -187,10 +198,10 @@ def create_boundary_condition_clipped_boundary( clip_box_args: tuple[Any, ...], ) -> Optional[StateType]: """ - Create a clipped boundary condition for a given state in the clipped model. + Create a boundary condition for the clipped model. + The function takes the original model as a reference to determine where - boundary conditions should NOT be placed, then applies this information to - create the boundary condition in the clipped model. + boundary conditions should NOT be placed. Parameters ---------- From fc1d283bd2c1033318ad88bf23e9d9e047688570 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 14:33:09 +0200 Subject: [PATCH 33/39] Rename cls to pkg_type here and add missing docstrings --- imod/mf6/utilities/imod5_converter.py | 107 ++++++++++++++++++++++++-- 1 file changed, 101 insertions(+), 6 deletions(-) diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index ef2f1d20f..0eac5e611 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -22,6 +22,22 @@ def convert_ibound_to_idomain( ibound: xr.DataArray, thickness: xr.DataArray ) -> xr.DataArray: + """ + Convert IBOUND array to IDOMAIN array. IBOUND -1 will be set to 1 in + IDOMAIN. When the thickness is <= 0, IDOMAIN will be set to -1. + + Parameters + ---------- + ibound : xr.DataArray + The IBOUND array from iMOD5. + thickness : xr.DataArray + The thickness array of the model layers. + + Returns + ------- + xr.DataArray + The corresponding IDOMAIN array. + """ # Convert IBOUND to IDOMAIN # -1 to 1, these will have to be filled with # CHD cells. @@ -127,6 +143,26 @@ def well_from_imod5_cap_data( ``"artificial_recharge_capacity"`` is ignored as the abstraction capacity is already defined in the point data. This is an ``n:1`` mapping: multiple grid cells can map to one well. + + Parameters + ---------- + imod5_data : Imod5DataDict + The iMOD5 data containing the "cap" package with abstraction + information. + target_dis : Optional[IRegridPackage] + The target discretization package for regridding the data. Required if + the data is in grid format (IDF). + regridder_types : DataclassType + The regrid methods to use for regridding the data. + regrid_cache : RegridderWeightsCache + Cache for storing regridder weights to speed up repeated regridding + operations. + + Returns + ------- + dict[str, np.ndarray] + A dictionary containing well information extracted from the iMOD5 cap + data. """ cap_data = imod5_data["cap"] has_ipf_well = isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame) @@ -146,7 +182,7 @@ def well_from_imod5_cap_data( def regrid_imod5_pkg_data( - cls: Optional[type[Package]], + pkg_type: Optional[type[Package]], imod5_pkg_data: GridDataDict, target_dis: Package, regridder_types: Optional[DataclassType], @@ -155,14 +191,37 @@ def regrid_imod5_pkg_data( """ Regrid iMOD5 package data to target idomain. Optionally get regrid methods from class if not provided. + + Parameters + ---------- + pkg_type: + The type of the package being regridded. This is used to determine the + appropriate regrid methods if regridder_types is not provided. + imod5_pkg_data: + The iMOD5 package data to be regridded. + target_dis: + The target discretization package containing the idomain to regrid to. + regridder_types: + Optional regrid methods to use for regridding. If not provided, they + will be obtained from the pkg_type. + regrid_cache: + Cache for storing regridder weights to speed up repeated regridding + operations. + + Returns + ------- + GridDataDict + The regridded iMOD5 package data. """ - if (cls is None) and (regridder_types is None): + if (pkg_type is None) and (regridder_types is None): raise ValueError( - "Either cls or regridder_types must be provided for regridding." + "Either pkg_type or regridder_types must be provided for regridding." ) # set up regridder methods - elif (cls is not None) and (regridder_types is None): # check cls not None for mypy - regridder_types = cls.get_regrid_methods() + elif (pkg_type is not None) and ( + regridder_types is None + ): # check pkg_type not None for mypy + regridder_types = pkg_type.get_regrid_methods() # For mypy to succeed regridder_types = cast(DataclassType, regridder_types) @@ -178,6 +237,22 @@ def regrid_imod5_pkg_data( def chd_cells_from_imod5_data( imod5_pkg_data: GridDataDict, target_idomain: GridDataArray ) -> GridDataDict: + """ + Get CHD cells from iMOD5 package data based on IBOUND and target idomain. + + Parameters + ---------- + imod5_pkg_data: + The iMOD5 package data containing "head" and "ibound". + target_idomain: + The target idomain to filter active cells. + + Returns + ------- + GridDataDict + The filtered CHD cells with "head" values where IBOUND < 0 and target + idomain > 0. + """ head = imod5_pkg_data["head"] ibound = imod5_pkg_data["ibound"] @@ -200,6 +275,26 @@ def mask_topsystem_packages_with_ibound( """ Mask all top system packages where IBOUND < 0. These locations are assigned a constant head. + + Parameters + ---------- + imod5_data: + The iMOD5 data containing the "bnd" package with "ibound". + model: + The target MODFLOW 6 model. + regridder_types: + Optional regrid methods to use for regridding. If not provided, default + methods will be used. + regrid_cache: + Cache for storing regridder weights to speed up repeated regridding + operations. + ignore_time_purge_empty: + Flag indicating whether to ignore time when purging empty cells. + + Returns + ------- + None + The function modifies the top system packages in the model in place. """ if regridder_types is None: @@ -207,7 +302,7 @@ def mask_topsystem_packages_with_ibound( ibound = imod5_data["bnd"]["ibound"] regridded_ibound = regrid_imod5_pkg_data( - cls=None, + pkg_type=None, imod5_pkg_data={"ibound": ibound}, target_dis=model["dis"], regridder_types=regridder_types, From 8061c3fe1c83c1789b7d8f22c40242e6e9f4a487 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 15:49:09 +0200 Subject: [PATCH 34/39] Move mask_topsyste_where_bc logic to utilities module --- imod/mf6/model.py | 14 +++----------- imod/mf6/utilities/mask.py | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 11 deletions(-) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index 6d783fd40..74cf37198 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -41,7 +41,7 @@ StateClassType, create_boundary_condition_clipped_boundary, ) -from imod.mf6.utilities.mask import mask_topsystem +from imod.mf6.utilities.mask import mask_topsystem_where_bc from imod.mf6.utilities.mf6hfb import merge_hfb_packages from imod.mf6.validation_settings import ValidationSettings from imod.mf6.wel import GridAgnosticWell @@ -730,17 +730,9 @@ def clip_box( clipped[pkg_name] = clipped_boundary_condition - # Mask topsystem packages where the state boundary cells have been - # added. - state_varname = clipped_boundary_condition._period_data[0] - # Select the state variable for the first time step as mask. - # Its location will be constant through time. - state_var = clipped_boundary_condition.dataset[state_varname].isel( - time=0, missing_dims="ignore", drop=True + mask_topsystem_where_bc( + clipped, clipped_boundary_condition, ignore_time_purge_empty ) - not_added_bc = np.isnan(state_var) - # Purge empty packages called by the mask_topsystem function - mask_topsystem(clipped, not_added_bc, ignore_time_purge_empty) else: clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty) diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py index d734d5518..5d0d1bd71 100644 --- a/imod/mf6/utilities/mask.py +++ b/imod/mf6/utilities/mask.py @@ -1,4 +1,9 @@ +from typing import cast + +import numpy as np + from imod.common.interfaces.imodel import IModel +from imod.common.interfaces.ipackage import IPackage from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition from imod.typing import GridDataArray @@ -26,3 +31,33 @@ def mask_topsystem( if isinstance(pkg, ITopSystemBoundaryCondition) ] model.mask_packages(topsystem_packages, is_active, ignore_time_purge_empty) + + +def mask_topsystem_where_bc( + model: IModel, boundary_condition: IPackage, ignore_time: bool +) -> None: + """ + Mask all top system packages in the model inplace where the boundary + condition is active. + + Parameters + ---------- + model : IModel + The MODFLOW 6 model containing top system packages. + boundary_condition : IPackage + The boundary condition package used to determine which cells are not + added. + ignore_time : bool + If True, ignore the time dimension when masking the packages and set the + mask to where the first time step of the boundary condition contains + active cells. Else, aggregate over all time steps to determine the mask. + """ + state_varname = cast(str, boundary_condition._period_data[0]) + state_var = boundary_condition.dataset[state_varname] + if "time" in state_var.dims: + if ignore_time: + state_var = state_var.isel(time=0, drop=True) + else: + state_var = state_var.min(dim="time") + not_added_bc = np.isnan(state_var) + mask_topsystem(model, not_added_bc, ignore_time) From a213c4617831664f23d23152b13ef6db0a37e2de Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 15:56:40 +0200 Subject: [PATCH 35/39] Add BoundaryCondition interface --- imod/common/interfaces/iboundarycondition.py | 12 ++++++++++++ imod/mf6/boundary_condition.py | 5 +++-- imod/mf6/utilities/mask.py | 10 ++++------ 3 files changed, 19 insertions(+), 8 deletions(-) create mode 100644 imod/common/interfaces/iboundarycondition.py diff --git a/imod/common/interfaces/iboundarycondition.py b/imod/common/interfaces/iboundarycondition.py new file mode 100644 index 000000000..7db2f9db2 --- /dev/null +++ b/imod/common/interfaces/iboundarycondition.py @@ -0,0 +1,12 @@ +from abc import ABC +from typing import ClassVar + +from imod.common.interfaces.ipackage import IPackage + + +class IBoundaryCondition(IPackage, ABC): + """ + Interface for boundary condition packages with stress-period variables. + """ + + _period_data: ClassVar[tuple[str, ...]] diff --git a/imod/mf6/boundary_condition.py b/imod/mf6/boundary_condition.py index 882ee10df..468c18dc6 100644 --- a/imod/mf6/boundary_condition.py +++ b/imod/mf6/boundary_condition.py @@ -7,6 +7,7 @@ import xarray as xr import xugrid as xu +from imod.common.interfaces.iboundarycondition import IBoundaryCondition from imod.common.utilities.value_filters import enforce_scalar from imod.mf6.auxiliary_variables import ( expand_transient_auxiliary_variables, @@ -57,7 +58,7 @@ def _disv_recarr(arrdict, layer, notnull): return recarr -class BoundaryCondition(Package, abc.ABC): +class BoundaryCondition(Package, IBoundaryCondition, abc.ABC): """ BoundaryCondition is used to share methods for specific stress packages with a time component. @@ -324,7 +325,7 @@ def _get_period_varnames(self) -> list[str]: >>> river._get_period_varnames() >>> # prints: ['stage', 'conductance', 'bottom_elevation', 'species1', 'species2'] """ - result = [] + result: list[str] = [] if hasattr(self, "_period_data"): result.extend(self._period_data) if hasattr(self, "_optional_data"): diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py index 5d0d1bd71..d2d6f13c8 100644 --- a/imod/mf6/utilities/mask.py +++ b/imod/mf6/utilities/mask.py @@ -1,9 +1,7 @@ -from typing import cast - import numpy as np +from imod.common.interfaces.iboundarycondition import IBoundaryCondition from imod.common.interfaces.imodel import IModel -from imod.common.interfaces.ipackage import IPackage from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition from imod.typing import GridDataArray @@ -34,7 +32,7 @@ def mask_topsystem( def mask_topsystem_where_bc( - model: IModel, boundary_condition: IPackage, ignore_time: bool + model: IModel, boundary_condition: IBoundaryCondition, ignore_time: bool ) -> None: """ Mask all top system packages in the model inplace where the boundary @@ -44,7 +42,7 @@ def mask_topsystem_where_bc( ---------- model : IModel The MODFLOW 6 model containing top system packages. - boundary_condition : IPackage + boundary_condition : IBoundaryCondition The boundary condition package used to determine which cells are not added. ignore_time : bool @@ -52,7 +50,7 @@ def mask_topsystem_where_bc( mask to where the first time step of the boundary condition contains active cells. Else, aggregate over all time steps to determine the mask. """ - state_varname = cast(str, boundary_condition._period_data[0]) + state_varname = boundary_condition._period_data[0] state_var = boundary_condition.dataset[state_varname] if "time" in state_var.dims: if ignore_time: From 1ea8cad369ef98e5690e98db6b11baa49cf35516 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 16:09:41 +0200 Subject: [PATCH 36/39] update forgotten mock patch --- imod/tests/test_mf6/test_mf6_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py index 42f19e85b..fc371b88c 100644 --- a/imod/tests/test_mf6/test_mf6_model.py +++ b/imod/tests/test_mf6/test_mf6_model.py @@ -255,7 +255,7 @@ def test_clip_box_without_state_for_boundary(self, model_type, pkg_type): pkg_id = pkg_type._pkg_id assert f"{pkg_id}_clipped" not in clipped - @mock.patch("imod.mf6.model.mask_topsystem") + @mock.patch("imod.mf6.model.mask_topsystem_where_bc") @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg") @pytest.mark.parametrize( "model_type, pkg_type", @@ -314,7 +314,7 @@ def test_clip_box_with_state_for_boundary( ) mask_topsystem_mock.assert_called_once() - @mock.patch("imod.mf6.model.mask_topsystem") + @mock.patch("imod.mf6.model.mask_topsystem_where_bc") @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg") @pytest.mark.parametrize( "model_type, pkg_type", From af386aefd8ccf512a735812258023e2a08e3c8a8 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 16:12:56 +0200 Subject: [PATCH 37/39] Clearer varnames in fixture --- imod/tests/fixtures/msw_imod5_cap_fixture.py | 52 ++++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py index 76bb895f6..3d46a3bda 100644 --- a/imod/tests/fixtures/msw_imod5_cap_fixture.py +++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py @@ -19,10 +19,10 @@ def imod5_cap_data() -> GridDataDict: da_kwargs["coords"] = {"layer": layer, "y": y, "x": x, "dx": dx, "dy": dy} imod5_data = {} - d = {} + d_cap = {} # fmt: off - d["boundary"] = xr.DataArray( + d_cap["boundary"] = xr.DataArray( np.array([ [ [1, 1, 1], @@ -32,7 +32,7 @@ def imod5_cap_data() -> GridDataDict: dtype=int), **da_kwargs ) - d["landuse"] = xr.DataArray( + d_cap["landuse"] = xr.DataArray( np.array([ [ [1, 2, 3], @@ -42,7 +42,7 @@ def imod5_cap_data() -> GridDataDict: dtype=int), **da_kwargs ) - d["rootzone_thickness"] = xr.DataArray( + d_cap["rootzone_thickness"] = xr.DataArray( np.array([ [ [0.1, 0.1, 0.1], @@ -52,7 +52,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["soil_physical_unit"] = xr.DataArray( + d_cap["soil_physical_unit"] = xr.DataArray( np.array([ [ [1, 1, 1], @@ -62,7 +62,7 @@ def imod5_cap_data() -> GridDataDict: dtype=int), **da_kwargs ) - d["surface_elevation"] = xr.DataArray( + d_cap["surface_elevation"] = xr.DataArray( np.array([ [ [1.1, 1.2, 1.3], @@ -72,7 +72,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["artificial_recharge"] = xr.DataArray( + d_cap["artificial_recharge"] = xr.DataArray( np.array([ [ [0.2, 0.2, 0.2], @@ -82,7 +82,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["artificial_recharge_layer"] = xr.DataArray( + d_cap["artificial_recharge_layer"] = xr.DataArray( np.array([ [ [1, 2, 3], @@ -92,7 +92,7 @@ def imod5_cap_data() -> GridDataDict: dtype=int), **da_kwargs ) - d["artificial_recharge_capacity"] = xr.DataArray( + d_cap["artificial_recharge_capacity"] = xr.DataArray( np.array([ [ [0.4, 0.4, 0.4], @@ -102,7 +102,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["wetted_area"] = xr.DataArray( + d_cap["wetted_area"] = xr.DataArray( np.array([ [ [0.1, 0.1, 0.1], @@ -112,7 +112,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["urban_area"] = xr.DataArray( + d_cap["urban_area"] = xr.DataArray( np.array([ [ [0.2, 0.2, 0.2], @@ -122,7 +122,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["urban_ponding_depth"] = xr.DataArray( + d_cap["urban_ponding_depth"] = xr.DataArray( np.array([ [ [2.2, 2.2, 2.2], @@ -132,7 +132,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["rural_ponding_depth"] = xr.DataArray( + d_cap["rural_ponding_depth"] = xr.DataArray( np.array([ [ [1.2, 1.2, 1.2], @@ -142,7 +142,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["urban_runoff_resistance"] = xr.DataArray( + d_cap["urban_runoff_resistance"] = xr.DataArray( np.array([ [ [3.2, 3.2, 3.2], @@ -152,7 +152,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["rural_runoff_resistance"] = xr.DataArray( + d_cap["rural_runoff_resistance"] = xr.DataArray( np.array([ [ [3.6, 3.6, 3.6], @@ -162,7 +162,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["urban_runon_resistance"] = xr.DataArray( + d_cap["urban_runon_resistance"] = xr.DataArray( np.array([ [ [5.2, 5.2, 5.2], @@ -172,7 +172,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["rural_runon_resistance"] = xr.DataArray( + d_cap["rural_runon_resistance"] = xr.DataArray( np.array([ [ [5.6, 5.6, 5.6], @@ -182,7 +182,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["urban_infiltration_capacity"] = xr.DataArray( + d_cap["urban_infiltration_capacity"] = xr.DataArray( np.array([ [ [10.2, 10.2, 10.2], @@ -192,7 +192,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["rural_infiltration_capacity"] = xr.DataArray( + d_cap["rural_infiltration_capacity"] = xr.DataArray( np.array([ [ [20.2, 20.2, 20.2], @@ -202,7 +202,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["perched_water_table_level"]= xr.DataArray( + d_cap["perched_water_table_level"]= xr.DataArray( np.array([ [ [2.0, 2.0, 2.0], @@ -212,7 +212,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["soil_moisture_fraction"]= xr.DataArray( + d_cap["soil_moisture_fraction"]= xr.DataArray( np.array([ [ [1.5, 1.5, 1.5], @@ -222,7 +222,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["conductivity_factor"]= xr.DataArray( + d_cap["conductivity_factor"]= xr.DataArray( np.array([ [ [2.5, 2.5, 2.5], @@ -232,8 +232,8 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d2 = {} - d2["ibound"] = xr.DataArray( + d_bnd = {} + d_bnd["ibound"] = xr.DataArray( np.array([ [ [1, 1, 1], @@ -244,7 +244,7 @@ def imod5_cap_data() -> GridDataDict: **da_kwargs ) # fmt: on - imod5_data["cap"] = d - imod5_data["bnd"] = d2 + imod5_data["cap"] = d_cap + imod5_data["bnd"] = d_bnd return imod5_data From dccb8b3406b0ec0f662ea2869f56e60bce8389cc Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 28 Sep 2026 16:14:00 +0200 Subject: [PATCH 38/39] Update imod/mf6/utilities/clipped_bc_creator.py Co-authored-by: Claire Donnelly <82878115+ClaireDons@users.noreply.github.com> --- imod/mf6/utilities/clipped_bc_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py index c666ae261..762f13e6f 100644 --- a/imod/mf6/utilities/clipped_bc_creator.py +++ b/imod/mf6/utilities/clipped_bc_creator.py @@ -245,7 +245,7 @@ def create_boundary_condition_clipped_boundary( # Remove all indices before first timestep of state_for_clipped_boundary. # This to prevent empty dataarrays unnecessarily being made for these # indices, which can lead to them to be removed when purging empty packages - # with ignore_time=True. Unfortunately, this is needs to be handled here and + # with ignore_time=True. Unfortunately, this needs to be handled here and # not in _create_boundary_condition_for_unassigned_boundary, as otherwise # this function is called twice which could result in broadcasting errors in # the second call if the time domain of state_for_boundary and assigned From ab7212fd27ecb0d8c2a02a1cd7de4382eb1562c1 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 28 Sep 2026 16:47:01 +0200 Subject: [PATCH 39/39] Add unittests --- .../test_utilities/test_mf6_mask_util.py | 70 ++++++++++++++++++- 1 file changed, 67 insertions(+), 3 deletions(-) diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py index 167d6d1d4..9edf9c549 100644 --- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py +++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py @@ -1,13 +1,16 @@ import numpy as np +import pandas as pd +import pytest_cases +import xarray as xr -from imod.mf6.utilities.mask import mask_topsystem +from imod.mf6.utilities.mask import mask_topsystem, mask_topsystem_where_bc from imod.typing.grid import ones_like, zeros_like def test_mask_topsystem(twri_model): """ - Test the mask_topsystem utility function by deactivating all cells in the - grid. + Test the mask_topsystem utility function by deactivating the first cell in + the grid. """ # Arrange gwf_model = twri_model["GWF_1"] @@ -38,3 +41,64 @@ def test_mask_topsystem__all_removed(twri_model): # Assert for key in ["rch", "drn"]: assert key not in gwf_model.keys() + + +def test_mask_topsystem_where_bc__no_time(twri_model): + """ + Test the mask_topsystem utility function with boundary conditions by deactivating the first cell in the grid. + """ + # Arrange + gwf_model = twri_model["GWF_1"] + bc_condition = gwf_model["chd"] + + # Test if fixture still as expected, first cell should be active in the + # boundary condition + chd_active = np.isnan(bc_condition.dataset["head"]) + assert chd_active[0, 0, 0].item() is False + assert chd_active[0, 0, 1].item() is True + # Should be all active + rch_active_prior = np.isnan(gwf_model["rch"].dataset["rate"]) + assert rch_active_prior[0, 0].item() is False + assert rch_active_prior[0, 1].item() is False + + # Act + mask_topsystem_where_bc(gwf_model, bc_condition, True) + + # Assert + rch_active_post = np.isnan(gwf_model["rch"].dataset["rate"]) + assert rch_active_post[0, 0].item() is True + assert rch_active_post[0, 1].item() is False + + +@pytest_cases.parametrize("ignore_time", [True, False]) +def test_mask_topsystem_where_bc__time(twri_model, ignore_time): + """ + Test the mask_topsystem utility function with boundary conditions by deactivating the first cell in the grid. + """ + # Arrange + gwf_model = twri_model["GWF_1"] + bc_condition = gwf_model["chd"] + + # Make package transient + time_coord = pd.date_range("2000-01-01", periods=3) + time_da = xr.DataArray([1, 1, 1], dims=["time"], coords={"time": time_coord}) + ds = bc_condition.dataset.copy() + bc_condition.dataset["head"] = time_da * ds["head"] + + # Test if fixture still as expected, first cell should be active in the + # boundary condition + chd_active = np.isnan(bc_condition.dataset["head"]) + assert chd_active[1, 0, 0, 0].item() is False + assert chd_active[1, 0, 0, 1].item() is True + # Should be all active + rch_active_prior = np.isnan(gwf_model["rch"].dataset["rate"]) + assert rch_active_prior[0, 0].item() is False + assert rch_active_prior[0, 1].item() is False + + # Act + mask_topsystem_where_bc(gwf_model, bc_condition, ignore_time) + + # Assert + rch_active_post = np.isnan(gwf_model["rch"].dataset["rate"]) + assert rch_active_post[0, 0].item() is True + assert rch_active_post[0, 1].item() is False