diff --git a/imod_coupler/drivers/ribametamod/exchange.py b/imod_coupler/drivers/ribametamod/exchange.py index 58389f93..c126660c 100644 --- a/imod_coupler/drivers/ribametamod/exchange.py +++ b/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( diff --git a/imod_coupler/drivers/ribametamod/mapping.py b/imod_coupler/drivers/ribametamod/mapping.py index 245ae740..86ae75ac 100644 --- a/imod_coupler/drivers/ribametamod/mapping.py +++ b/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: diff --git a/imod_coupler/drivers/ribametamod/ribametamod.py b/imod_coupler/drivers/ribametamod/ribametamod.py index c7bbb29d..3abf4fdc 100644 --- a/imod_coupler/drivers/ribametamod/ribametamod.py +++ b/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, }, ) ) diff --git a/tests/fixtures/fixture_ribasim.py b/tests/fixtures/fixture_ribasim.py index add7bf4b..0bc576ba 100644 --- a/tests/fixtures/fixture_ribasim.py +++ b/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()) diff --git a/tests/test_imod_coupler/test_ribametamod.py b/tests/test_imod_coupler/test_ribametamod.py index 653fc77a..19029185 100644 --- a/tests/test_imod_coupler/test_ribametamod.py +++ b/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( diff --git a/tests/test_imod_coupler/test_ribametamod_cases.py b/tests/test_imod_coupler/test_ribametamod_cases.py index a037859e..39dda5c2 100644 --- a/tests/test_imod_coupler/test_ribametamod_cases.py +++ b/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,