From 200fb59224b7735fb215baf70260a3dd5c4a944a Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Mon, 5 Feb 2024 15:27:35 +0100 Subject: [PATCH 1/3] Add check for time window on models --- pre-processing/primod/ribamod/ribamod.py | 33 ++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/pre-processing/primod/ribamod/ribamod.py b/pre-processing/primod/ribamod/ribamod.py index 1c9971f3..b57d09b0 100644 --- a/pre-processing/primod/ribamod/ribamod.py +++ b/pre-processing/primod/ribamod/ribamod.py @@ -8,6 +8,7 @@ import pandas as pd import ribasim import tomli_w +import xarray as xr from imod.mf6 import Drainage, GroundwaterFlowModel, Modflow6Simulation, River from primod.ribamod.exchange_creator import ( @@ -65,6 +66,10 @@ def __init__( coupling_list: list[DriverCoupling], basin_definition: gpd.GeoDataFrame, ): + self.validate_time_window( + ribasim_model=ribasim_model, + mf6_simulation=mf6_simulation, + ) self.ribasim_model = ribasim_model self.mf6_simulation = mf6_simulation self.coupling_list = coupling_list @@ -202,6 +207,33 @@ def validate_keys( f"present in the model: {missing}" ) + @staticmethod + def validate_time_window( + ribasim_model: ribasim.Model, + mf6_simulation: Modflow6Simulation, + ) -> None: + def to_timestamp(xr_time: xr.DataArray) -> pd.Timestamp: + return pd.Timestamp(xr_time.to_numpy().item()) + + mf6_timedis = mf6_simulation["time_discretization"].dataset + mf6_start = to_timestamp(mf6_timedis["time"].isel(time=0)).to_pydatetime() + time_delta = pd.to_timedelta( + mf6_timedis["timestep_duration"].isel(time=-1).item(), unit="days" + ) + mf6_end = ( + to_timestamp(mf6_timedis["time"].isel(time=-1)) + time_delta + ).to_pydatetime() + + ribasim_start = ribasim_model.starttime + ribasim_end = ribasim_model.endtime + if ribasim_start != mf6_start or ribasim_end != mf6_end: + raise ValueError( + "Ribasim simulation time window does not match MODFLOW6.\n" + f"Ribasim: {ribasim_start} to {ribasim_end}\n" + f"MODFLOW6: {mf6_start} to {mf6_end}\n" + ) + return + def validate_basin_node_ids(self) -> pd.Series: assert self.ribasim_model.basin.profile.df is not None basin_ids = np.unique(self.ribasim_model.basin.profile.df["node_id"]) @@ -234,6 +266,7 @@ def write_exchanges( # Assume only one groundwater flow model # FUTURE: Support multiple groundwater flow models. + # FUTURE: move validations to the init? gwf_model = self.mf6_simulation[gwf_names[0]] self.validate_keys( gwf_model, From c724a6cedd3c82e78ced5b727a153c932598cb3e Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Mon, 5 Feb 2024 16:04:19 +0100 Subject: [PATCH 2/3] harmonize dates with ribasim-testmodels --- tests/fixtures/fixture_modflow.py | 7 ++++++- tests/test_primod/test_preprocessing_ribamod.py | 10 ++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/tests/fixtures/fixture_modflow.py b/tests/fixtures/fixture_modflow.py index e8aa05ea..374ace87 100644 --- a/tests/fixtures/fixture_modflow.py +++ b/tests/fixtures/fixture_modflow.py @@ -176,6 +176,11 @@ def mf6_bucket_model(active_idomain: xr.DataArray) -> mf6.Modflow6Simulation: ) simulation = make_mf6_simulation(gwf_model) + # Override time discretization: + # Make model returns dates in 1970 + # Ribasim bucket models is 2020-2021. + times = pd.date_range("2020-01-01", "2021-01-01", freq="D") + simulation.create_time_discretization(additional_times=times) return simulation @@ -286,7 +291,7 @@ def mf6_backwater_model() -> mf6.Modflow6Simulation: reordering_method=None, relaxation_factor=0.97, ) - times = pd.date_range("2020-01-01", "2030-01-01", freq="M") + times = pd.date_range("2020-01-01", "2021-01-01", freq="D") simulation.create_time_discretization(additional_times=times) return simulation diff --git a/tests/test_primod/test_preprocessing_ribamod.py b/tests/test_primod/test_preprocessing_ribamod.py index 7b306a90..cee36006 100644 --- a/tests/test_primod/test_preprocessing_ribamod.py +++ b/tests/test_primod/test_preprocessing_ribamod.py @@ -37,6 +37,16 @@ def basin_definition(mf6_bucket_model, ribasim_bucket_model) -> gpd.GeoDataFrame return gpd.GeoDataFrame(data={"node_id": node_id}, geometry=[polygon]) +def test_validate_time_window(mf6_bucket_model, ribasim_bucket_model): + times = pd.date_range("1970-01-01", "1971-01-01", freq="D") + # override time discretization + mf6_bucket_model.create_time_discretization(additional_times=times) + with pytest.raises( + ValueError, match="Ribasim simulation time window does not match MODFLOW6" + ): + RibaMod.validate_time_window(ribasim_bucket_model, mf6_bucket_model) + + def test_validate_keys(mf6_bucket_model): with pytest.raises(ValueError, match="active and passive keys share members"): RibaMod.validate_keys( From 2fb0a752807a86f3079ff8295056321d05bf0551 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Mon, 5 Feb 2024 20:51:33 +0100 Subject: [PATCH 3/3] Adjust tests for updated testmodels (shorter simulation periods) --- tests/fixtures/fixture_modflow.py | 2 +- tests/test_ribamod.py | 15 +++++---------- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/tests/fixtures/fixture_modflow.py b/tests/fixtures/fixture_modflow.py index 374ace87..6a06e307 100644 --- a/tests/fixtures/fixture_modflow.py +++ b/tests/fixtures/fixture_modflow.py @@ -333,7 +333,7 @@ def two_basin_variation( reordering_method=None, relaxation_factor=0.97, ) - times = pd.date_range("2020-01-01", "2020-02-01", freq="1d") + times = pd.date_range("2020-01-01", "2021-01-01", freq="1d") simulation.create_time_discretization(additional_times=times) return simulation diff --git a/tests/test_ribamod.py b/tests/test_ribamod.py index a5a26291..20407cc3 100644 --- a/tests/test_ribamod.py +++ b/tests/test_ribamod.py @@ -7,7 +7,6 @@ import pandas as pd import pytest import xarray as xr -from numpy.testing import assert_allclose from primod.ribamod import RibaMod from pytest_cases import parametrize_with_cases @@ -136,22 +135,18 @@ def test_ribamod_backwater( run_coupler_function=run_coupler_function, ) - final_level = results.basin_df[results.basin_df["time"] == "2029-12-01"]["level"] + final_level = results.basin_df[results.basin_df["time"] == "2021-01-01"]["level"] # Assert that the final level is a mototonically decreasing curve. assert (np.diff(final_level) < 0).all() # The head should follow the same pattern. assert (results.head.isel(layer=0, time=-1).diff("x") < 0.0).all() - drn = results.budgets["drn-1"].compute() - riv = results.budgets["riv-1"].compute() - # At the last time step, the drain and the river should have equal water - # balance terms since they have the same conductance and the river_stage = - # drainage_elevation. - assert_allclose(drn.isel(time=-1), riv.isel(time=-1)) + drn = results.budgets["drn-1"].isel(time=-1).compute() + riv = results.budgets["riv-1"].isel(time=-1).compute() # Get the last flow between the edges - final_flow = results.flow_df[results.flow_df["time"] == "2029-12-01"] + final_flow = results.flow_df[results.flow_df["time"] == "2021-01-01"] # Check's what lost and gained in the basins network = ribamod_model.ribasim_model.network basin_ids = network.node.df.index[network.node.df["type"] == "Basin"] @@ -162,7 +157,7 @@ def test_ribamod_backwater( ]["flow"] * 86_400.0 ).to_numpy() - modflow_budget = (drn + riv).isel(time=-1).sel(y=0).to_numpy() + modflow_budget = (drn + riv).sel(y=0).to_numpy() budget_diff = ribasim_budget + modflow_budget assert (np.abs(budget_diff) < 1.0e-6).all()