Test ribasim integration (#334)

Improvement of the volume integration for realised infiltration and drainage fluxes on the Ribasim side

.teamcity/_Self/buildTypes/TestbenchCouplerWin64_2.kt

@@ -92,6 +92,10 @@ object TestbenchCouplerWin64_2 : BuildType({
         }
     }
 
+    failureConditions {
+        executionTimeoutMin = 120
+    }
+
     dependencies {
         dependency(IMODCollector.buildTypes.IMODCollector_X64development) {
             snapshot {

imod_coupler/drivers/ribametamod/ribametamod.py

@@ -73,9 +73,12 @@ class RibaMetaMod(Driver):
     ribasim_level: NDArray[Any]
     ribasim_infiltration: NDArray[Any]
     ribasim_drainage: NDArray[Any]
+    ribasim_infiltration_save: NDArray[Any]
+    ribasim_drainage_save: NDArray[Any]
     ribasim_volume: NDArray[Any]
     ribasim_user_demand: NDArray[Any]
     ribasim_user_realized: NDArray[Any]
+    ribasim_user_realized_save: NDArray[Any]
 
     # MetaSWAP variables
     mf6_sprinkling_wells: NDArray[Any]  # the well data for coupled extractions
@@ -197,18 +200,23 @@ def couple_ribasim(self, mf6_flowmodel_key: str) -> ChainMap[str, Any]:
             # Get all Ribasim pointers, relevant for coupling with MODFLOW 6
             self.ribasim_infiltration = self.ribasim.get_value_ptr("basin.infiltration")
             self.ribasim_drainage = self.ribasim.get_value_ptr("basin.drainage")
-            self.ribasim_infiltration_sum = self.ribasim.get_value_ptr(
+            self.ribasim_cumulative_infiltration = self.ribasim.get_value_ptr(
                 "basin.cumulative_infiltration"
             )
-            self.ribasim_drainage_sum = self.ribasim.get_value_ptr(
+            self.ribasim_infiltration_save = np.empty_like(
+                self.ribasim_cumulative_infiltration
+            )
+            self.ribasim_cumulative_drainage = self.ribasim.get_value_ptr(
                 "basin.cumulative_drainage"
             )
+            self.ribasim_drainage_save = np.empty_like(self.ribasim_cumulative_drainage)
             self.ribasim_level = self.ribasim.get_value_ptr("basin.level")
             self.ribasim_storage = self.ribasim.get_value_ptr("basin.storage")
             self.ribasim_user_demand = self.ribasim.get_value_ptr("user_demand.demand")
             self.ribasim_user_realized = self.ribasim.get_value_ptr(
-                "user_demand.inflow"
+                "user_demand.cumulative_inflow"
             )
+            self.ribasim_user_realized_save = np.empty_like(self.ribasim_user_realized)
             self.subgrid_level = self.ribasim.get_value_ptr("basin.subgrid_level")
 
             # add to return ChainMap
@@ -364,9 +372,6 @@ def update_ribasim_metaswap(self) -> None:
             self.exchange.add_ponding_volume_msw(self.msw_ponding_volume)
             if self.enable_sprinkling_surface_water:
                 self.exchange_sprinkling_demand_msw2rib()
-                self.ribasim_user_realized[:] = (
-                    0.0  # reset cummulative for the next timestep
-                )
             # exchange summed volumes to Ribasim
             self.exchange.flux_to_ribasim(self.delt_gw, self.delt_sw)
             # update Ribasim per delt_sw
@@ -399,8 +404,13 @@ def update(self) -> None:
                 self.update_ribasim_metaswap()
             else:
                 self.update_ribasim()
+
             self.exchange.flux_to_modflow(
-                (self.ribasim_drainage_sum - self.ribasim_infiltration_sum),
+                (self.ribasim_cumulative_drainage[:] - self.ribasim_drainage_save[:])
+                - (
+                    self.ribasim_cumulative_infiltration[:]
+                    - self.ribasim_infiltration_save[:]
+                ),
                 self.delt_gw,
             )
             self.exchange.log_demands(self.get_current_time())
@@ -464,8 +474,8 @@ def exchange_rib2mod(self) -> None:
     def exchange_mod2rib(self) -> None:
         self.exchange.add_flux_estimate_mod(self.mf6_head, self.delt_gw)
         # reset Ribasim pointers
-        self.ribasim_infiltration_sum[:] = 0.0
-        self.ribasim_drainage_sum[:] = 0.0
+        self.ribasim_infiltration_save[:] = self.ribasim_cumulative_infiltration[:]
+        self.ribasim_drainage_save[:] = self.ribasim_cumulative_drainage[:]
 
     def exchange_sprinkling_demand_msw2rib(self) -> None:
         # flux demand from metaswap sprinkling to Ribasim (demand)
@@ -494,7 +504,10 @@ def exchange_sprinkling_flux_realised_msw2rib(self) -> None:
 
         self.realised_fractions_swspr[:] = 1.0  # all realized for non-coupled svats
         self.realised_fractions_swspr[nonzero_user_indices] = (
-            self.ribasim_user_realized[nonzero_user_indices]
+            (
+                self.ribasim_user_realized[nonzero_user_indices]
+                - self.ribasim_user_realized_save[nonzero_user_indices]
+            )
             / (self.delt_sw * day_to_seconds)
         ) / self.mapped_sprinkling_demand[nonzero_user_indices]
 
@@ -504,7 +517,9 @@ def exchange_sprinkling_flux_realised_msw2rib(self) -> None:
         msw_sprinkling_realized[:] = (
             self.msw.get_surfacewater_sprinking_demand_ptr() * msw_sprfrac_realised
         )[:]
-        self.ribasim_user_realized[:] = 0.0  # reset cummulative for the next timestep
+        self.ribasim_user_realized_save[:] = self.ribasim_user_realized[
+            :
+        ]  # keep old values
         self.exchange_logger.log_exchange(
             ("sprinkling_realized"),
             msw_sprinkling_realized,

tests/fixtures/fixture_ribasim.py

@@ -6,6 +6,8 @@
 import ribasim
 import ribasim_testmodels
 
+solver_algorithm: str = "QNDF"
+
 
 def add_subgrid(model: ribasim.Model) -> ribasim.Model:
     """Add 1:1 subgrid levels to model"""
@@ -33,29 +35,36 @@ def add_subgrid(model: ribasim.Model) -> ribasim.Model:
 def ribasim_bucket_model() -> ribasim.Model:
     bucket = ribasim_testmodels.bucket_model()
     bucket.endtime = datetime(2023, 1, 1, 0, 0)
+    bucket.solver.algorithm = solver_algorithm
     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)
+    bucket.solver.algorithm = solver_algorithm
     return bucket
 
 
 @pytest_cases.fixture(scope="function")
 def ribasim_backwater_model() -> ribasim.Model:
-    return add_subgrid(ribasim_testmodels.backwater_model())
+    backwater = ribasim_testmodels.backwater_model()
+    backwater.solver.algorithm = solver_algorithm
+    return add_subgrid(backwater)
 
 
 @pytest_cases.fixture(scope="function")
 def ribasim_two_basin_model() -> ribasim.Model:
-    return ribasim_testmodels.two_basin_model()
+    twobasin = ribasim_testmodels.two_basin_model()
+    twobasin.solver.algorithm = solver_algorithm
+    return twobasin
 
 
 @pytest_cases.fixture(scope="function")
 def ribasim_two_basin_model_dbg() -> ribasim.Model:
     model = ribasim_testmodels.two_basin_model()
     #   model.logging.verbosity = ribasim.Verbosity("debug")
     #   model.logging.verbosity = "debug"
+    model.solver.algorithm = solver_algorithm
     return model