fix issue and add test (#328)

subgrid for ribasim optional

imod_coupler/drivers/ribametamod/exchange.py

@@ -187,17 +187,18 @@ def flux_to_ribasim(self, delt_gw: float, delt_sw: float) -> None:
  demand_per_subtimestep = self.get_demand_flux_sec(delt_gw, delt_sw)
 
  # exchange to Ribasim; negative demand in exchange class means infiltration from Ribasim
- coupled_index = self.mapping.coupled_mod2rib
- self.ribasim_infiltration[coupled_index] = np.where(
+ self.ribasim_infiltration[self.mapping.coupled_index] = np.where(
  demand_per_subtimestep < 0, -demand_per_subtimestep, 0
- )[coupled_index]
- self.ribasim_drainage[coupled_index] = np.where(
+ )[self.mapping.coupled_index]
+ self.ribasim_drainage[self.mapping.coupled_index] = np.where(
  demand_per_subtimestep > 0, demand_per_subtimestep, 0
- )[coupled_index]
+ )[self.mapping.coupled_index]
 
  def flux_to_modflow(
  self, realised_volume: NDArray[np.float64], delt_gw: float
  ) -> None:
+ if not self.mf6_river_packages:
+ return # no active coupling
  super().compute_realised(realised_volume)
  for key in self.mf6_active_river_api_packages.keys():
  realised_fraction = np.where(

imod_coupler/drivers/ribametamod/mapping.py

@@ -32,10 +32,12 @@ def __init__(
  self.coupling = coupling
  if has_ribasim:
  self.set_ribasim_modflow_mapping(packages)
+ self.coupled_index = self.coupled_mod2rib
  if has_metaswap and mod2svat is not None:
  self.set_metaswap_modflow_mapping(packages, mod2svat)
  if has_ribasim and has_metaswap and mod2svat is not None:
  self.set_metaswap_ribasim_mapping(packages, mod2svat)
+ self.coupled_index |= self.coupled_msw2rib
 
  def set_ribasim_modflow_mapping(self, packages: ChainMap[str, Any]) -> None:
  self.map_mod2rib = {}
@@ -200,6 +202,9 @@ def set_metaswap_ribasim_mapping(
  self, packages: ChainMap[str, Any], mod2svat: Path
  ) -> None:
  table_node2svat: NDArray[np.int32]
+ self.coupled_msw2rib: NDArray[np.bool_] = np.full(
+ packages["ribasim_nbasin"], False
+ )
  self.msw2rib = {}
 
  # surface water ponding mapping
@@ -222,6 +227,7 @@ def set_metaswap_ribasim_mapping(
  packages["ribasim_nbasin"],
  "sum",
  )
+ self.coupled_msw2rib |= self.msw2rib["sw_ponding"].getnnz(axis=1) > 0
 
  # surface water sprinkling mapping
  if self.coupling.rib_msw_sprinkling_map_surface_water is not None:

imod_coupler/drivers/ribametamod/ribametamod.py

@@ -218,12 +218,12 @@ def couple_ribasim(self, mf6_flowmodel_key: str) -> ChainMap[str, Any]:
  self.mf6_drainage_packages,
  {
  "ribasim_nbasin": len(self.ribasim_level),
- "ribasim_nuser": (
-  len(self.ribasim_user_realized)
-  if self.ribasim_user_realized.ndim > 0
-  else 0
- ),
- "ribasim_nsubgrid": len(self.subgrid_level),
+ "ribasim_nuser": len(self.ribasim_user_realized)
+ if self.ribasim_user_realized.ndim > 0
+ else 0,
+ "ribasim_nsubgrid": len(self.subgrid_level)
+ if self.subgrid_level.ndim > 0
+ else 0,
  },
  )
  )

tests/fixtures/fixture_ribasim.py

@@ -36,6 +36,13 @@ def ribasim_bucket_model() -> ribasim.Model:
  return add_subgrid(bucket)
 
 
+@pytest_cases.fixture(scope="function")
+def ribasim_bucket_model_no_subgrid() -> ribasim.Model:
+ bucket = ribasim_testmodels.bucket_model()
+ bucket.endtime = datetime(2023, 1, 1, 0, 0)
+ return bucket
+
+
 @pytest_cases.fixture(scope="function")
 def ribasim_backwater_model() -> ribasim.Model:
  return add_subgrid(ribasim_testmodels.backwater_model())

tests/test_imod_coupler/test_ribametamod.py

@@ -545,7 +545,6 @@ def test_ribametamod_backwater(
 
 @pytest.mark.xdist_group(name="ribasim")
 @parametrize_with_cases("ribametamod_model", glob="bucket_model")
-# @pytest.mark.skip("potentially hangs, to be solved in model")
 def test_ribametamod_bucket(
  tmp_path_dev: Path,
  ribametamod_model: RibaMetaMod,
@@ -577,6 +576,34 @@ def test_ribametamod_bucket(
  assert_results(tmp_path_dev, ribametamod_model, results)
 
 
+@pytest.mark.xdist_group(name="ribasim")
+@parametrize_with_cases("ribametamod_model", glob="bucket_model_no_subgrid")
+def test_ribametamod_bucket_no_subgrid(
+ tmp_path_dev: Path,
+ ribametamod_model: RibaMetaMod,
+ modflow_dll_devel: Path,
+ ribasim_dll_devel: Path,
+ ribasim_dll_dep_dir_devel: Path,
+ metaswap_dll_devel: Path,
+ metaswap_dll_dep_dir_devel: Path,
+ run_coupler_function: Callable[[Path], None],
+) -> None:
+ """
+ Test if the bucket model runs without a subgrid in the Ribasim model
+ """
+ ribametamod_model.write(
+ tmp_path_dev,
+ modflow6_dll=modflow_dll_devel,
+ ribasim_dll=ribasim_dll_devel,
+ ribasim_dll_dependency=ribasim_dll_dep_dir_devel,
+ metaswap_dll=metaswap_dll_devel,
+ metaswap_dll_dependency=metaswap_dll_dep_dir_devel,
+ modflow6_write_kwargs={"binary": False},
+ )
+
+ run_coupler_function(tmp_path_dev / ribametamod_model._toml_name)
+
+
 @pytest.mark.xdist_group(name="ribasim")
 @parametrize_with_cases("ribametamod_model", glob="two_basin_model")
 def test_ribametamod_two_basin(

tests/test_imod_coupler/test_ribametamod_cases.py

@@ -7,6 +7,7 @@
 from fixtures.common import create_wells_max_layer
 from imod import msw
 from imod.mf6 import (
+ Drainage,
  Modflow6Simulation,
  Recharge,
  StorageCoefficient,
@@ -17,6 +18,7 @@
  RibaMetaDriverCoupling,
  RibaMetaMod,
  RibaModActiveDriverCoupling,
+ RibaModPassiveDriverCoupling,
 )
 from ribasim.geometry import NodeTable as RibaNodeTbl
 from ribasim.nodes import level_demand, user_demand
@@ -326,6 +328,59 @@ def case_bucket_model(
  )
 
 
+def case_bucket_model_no_subgrid(
+ mf6_bucket_model: Modflow6Simulation,
+ msw_bucket_model: MetaSwapModel,
+ ribasim_bucket_model_no_subgrid: ribasim.Model,
+) -> RibaMetaMod:
+ mf6_bucket_model = set_confined_storage_formulation(mf6_bucket_model)
+
+ # add rch and well-package for coupling MetaMod
+ mf6_bucket_model = add_rch_package(mf6_bucket_model)
+ mf6_bucket_model = add_well_package(mf6_bucket_model)
+
+ # get variables
+ mf6_modelname, mf6_model = get_mf6_gwf_modelnames(mf6_bucket_model)[0]
+ mf6_active_river_packages = get_mf6_river_packagenames(mf6_model)
+ basin_definition = create_basin_definition(
+ ribasim_bucket_model_no_subgrid.basin.node, buffersize=100.0
+ )
+
+ # riv to drn package
+ conductance = (
+ mf6_bucket_model["GWF_1"][mf6_active_river_packages[0]].dataset["conductance"]
+ / 400
+ )
+ stage = mf6_bucket_model["GWF_1"][mf6_active_river_packages[0]].dataset[
+ "conductance"
+ ]
+ mf6_bucket_model["GWF_1"].pop(mf6_active_river_packages[0])
+ mf6_bucket_model["GWF_1"][mf6_active_river_packages[0]] = Drainage(
+ elevation=stage, conductance=conductance
+ )
+
+ metamod_coupling = MetaModDriverCoupling(
+ mf6_model=mf6_modelname,
+ mf6_recharge_package="rch_msw",
+ mf6_wel_package="well_msw",
+ )
+ ribamod_coupling = RibaModPassiveDriverCoupling(
+ mf6_model=mf6_modelname,
+ mf6_packages=mf6_active_river_packages,
+ ribasim_basin_definition=basin_definition,
+ )
+ ribameta_coupling = RibaMetaDriverCoupling(
+ ribasim_basin_definition=basin_definition,
+ )
+
+ return RibaMetaMod(
+ ribasim_model=ribasim_bucket_model_no_subgrid,
+ msw_model=msw_bucket_model,
+ mf6_simulation=mf6_bucket_model,
+ coupling_list=[metamod_coupling, ribamod_coupling, ribameta_coupling],
+ )
+
+
 def case_backwater_model(
  mf6_backwater_model: Modflow6Simulation,
  msw_backwater_model: MetaSwapModel,