Lumbricus update test-f + fix issues

imod_coupler/drivers/dfm_metamod/dfm_metamod.py

@@ -189,7 +189,7 @@ def update(self) -> None:
             self.exchange_stage_2d_dfm2msw()
 
             # initiate surface water timestep
-            self.msw.start_surface_water_time_step(idtsw)
+            self.msw.start_surface_water_time_step(idtsw + 1)
 
             # flux from metaswap ponding to water balance 1d
             self.exchange_ponding_msw2dflow1d()
@@ -219,11 +219,7 @@ def update(self) -> None:
             )
 
             # run dflow
-            while (
-                self.dfm.get_current_time()
-                < days_to_seconds(subtimestep_endtime) - self.time_eps
-            ):
-                self.dfm.update()
+            self.dfm.update(dt=days_to_seconds(self.delt_msw_dflow))
 
             # get cummelative flux after dfm-run
             time_after = self.dfm.get_current_time()
@@ -257,12 +253,10 @@ def update(self) -> None:
             # exchange 2d stage to msw, so it can finish the sw-timestep (now stage at the start of next timestep)
             self.exchange_stage_2d_dfm2msw()
 
-            self.msw.finish_surface_water_time_step(idtsw)
+            self.msw.finish_surface_water_time_step(idtsw + 1)
 
         # exchange correction flux to MF6
-        self.exchange_correction_dflow2mf6(
-            self.exchange_balans_1d.realised["dflow1d_flux2mf-riv_negative"]
-        )
+        self.exchange_correction_dflow2mf6()
 
         # convergence loop modflow-metaswap
         self.mf6.prepare_solve(1)
@@ -301,7 +295,7 @@ def get_array_dims(self) -> None:
             "mf6_recharge": self.mf6.get_recharge(
                 self.coupling.mf6_model, self.coupling.mf6_msw_recharge_pkg
             ).size,
-            "mf6_riv_active": self.mf6.get_river_flux_estimate(
+            "mf6_riv_active": self.mf6.get_river_drain_flux_estimate(
                 self.coupling.mf6_model, self.coupling.mf6_river_active_pkg
             ).size,
             "mf6_riv_passive": self.mf6.get_river_drain_flux(
@@ -517,7 +511,7 @@ def exchange_flux_riv_active_mf62dfm(self) -> None:
         """
         if self.map_active_mod_dflow1d["mf-riv2dflow1d_flux"] is not None:
             # conversion from (-)m3/dtgw to (+)m3/s
-            self.mf6_river_aquifer_flux_day = self.mf6.get_river_flux_estimate(
+            self.mf6_river_aquifer_flux_day = self.mf6.get_river_drain_flux_estimate(
                 self.coupling.mf6_model, self.coupling.mf6_river_active_pkg
             )
             mf6_river_aquifer_flux_sec = (
@@ -634,9 +628,14 @@ def exchange_flux_riv_passive_mf62dfm(self) -> None:
         MF6 unit: m3/d
         DFM unit: m3/s
         """
-        mf6_riv2_flux = self.mf6.get_river_drain_flux(
+        # mf6_riv2_flux = self.mf6.get_river_drain_flux(
+        #     self.coupling.mf6_model, self.coupling.mf6_river_passive_pkg
+        # )
+
+        mf6_riv2_flux = self.mf6.get_river_drain_flux_estimate(
             self.coupling.mf6_model, self.coupling.mf6_river_passive_pkg
         )
+
         # conversion from (-)m3/dtgw to (+)m3/s
         mf6_riv2_flux_sec = -mf6_riv2_flux / days_to_seconds(self.delt_mf6)
 
@@ -666,9 +665,14 @@ def exchange_flux_drn_passive_mf62dfm(
         MF6 unit: m3/d
         DFM unit: m3/s
         """
-        mf6_drn_flux = self.mf6.get_river_drain_flux(
-            self.coupling.mf6_model, self.coupling.mf6_drain_pkg
+        # mf6_drn_flux = self.mf6.get_river_drain_flux(
+        #     self.coupling.mf6_model, self.coupling.mf6_drain_pkg
+        # )
+
+        mf6_drn_flux = self.mf6.get_river_drain_flux_estimate(
+            self.coupling.mf6_model, self.coupling.mf6_drain_pkg, isdrain=True
         )
+
         # conversion from (+)m3/dtgw to (+)m3/s
         mf6_drn_flux_sec = -mf6_drn_flux / days_to_seconds(self.delt_mf6)
 
@@ -686,7 +690,7 @@ def exchange_flux_drn_passive_mf62dfm(
             self.dfm.get_current_time_days(),
         )
 
-    def exchange_correction_dflow2mf6(self, qdfm_realised: NDArray[float_]) -> None:
+    def exchange_correction_dflow2mf6(self) -> None:
         """
         From DFM to MF6
         the drainage/inflitration flux to the 1d rivers as realised by DFM is passed to
@@ -695,14 +699,11 @@ def exchange_correction_dflow2mf6(self, qdfm_realised: NDArray[float_]) -> None:
 
         if self.map_active_mod_dflow1d["mf-riv2dflow1d_flux"] is not None:
             wbal = self.exchange_balans_1d
-            realised_dfm = wbal.realised["dflow1d_flux2mf-riv_negative"]
-            demand_pos = wbal.demand["mf-riv2dflow1d_flux_positive"]
+            realised_neg = wbal.realised["dflow1d_flux2mf-riv_negative"]
             demand_neg = wbal.demand["mf-riv2dflow1d_flux_negative"]
             mask = np.nonzero(demand_neg)  # prevent zero division
-            realised_fraction = realised_dfm * 0.0 + 1.0
-            realised_fraction[mask] = (
-                realised_dfm[mask] - demand_pos[mask]
-            ) / demand_neg[mask]
+            realised_fraction = realised_neg * 0.0 + 1.0
+            realised_fraction[mask] = realised_neg[mask] / demand_neg[mask]
             matrix = self.map_active_mod_dflow1d["mf-riv2dflow1d_flux"].transpose()
 
             # correction only applies to Modflow cells which negatively contribute to the dflowfm volumes

imod_coupler/drivers/dfm_metamod/exchange.py

@@ -217,12 +217,15 @@ def check_maximum_shortage(
         sum_to_dflow : np.float_
             flux send to dflow
         """
-
-        shortage = np.absolute(sum_to_dflow - sum_from_dflow)
-        demand_msw = self.demand["msw-sprinkling2dflow1d_flux"]
-        demand_mf6_negative = self.demand["mf-riv2dflow1d_flux_negative"]
+        id = 10
+        shortage = np.around(np.absolute(sum_to_dflow - sum_from_dflow), decimals=id)
+        demand = np.around(
+            self.demand["msw-sprinkling2dflow1d_flux"]
+            + self.demand["mf-riv2dflow1d_flux_negative"],
+            decimals=id,
+        )
         condition = np.logical_and(
-            shortage > np.absolute(demand_msw + demand_mf6_negative),
+            shortage > np.absolute(demand),
             np.less(sum_from_dflow.astype(np.float32), sum_to_dflow.astype(np.float32)),
         )
         if np.any(condition):
@@ -240,7 +243,15 @@ def compute_realised(self, sum_from_dflow: NDArray[np.float_]) -> None:
            array with summed realised fluxes per dtsw-timstep by dflow
         """
 
+        # initialize realised arrays for cases where all conditions are false due to rounding errors
+        self.realised["dflow1d_flux2sprinkling_msw"] = self.demand[
+            "msw-sprinkling2dflow1d_flux"
+        ].copy()
+        self.realised["dflow1d_flux2mf-riv_negative"] = self.demand[
+            "mf-riv2dflow1d_flux_negative"
+        ].copy()
         sum_to_dflow = self.demand["sum"][:]
+
         # update elements for no shortage
         self.set_realised_no_shortage(sum_from_dflow, sum_to_dflow)
         # update elements for cases with shortage + shortage <= msw_demand
@@ -249,3 +260,5 @@ def compute_realised(self, sum_from_dflow: NDArray[np.float_]) -> None:
         self.set_realised_shortage_msw_mf6(sum_from_dflow, sum_to_dflow)
         # check if shortage is not larger than negative demand
         self.check_maximum_shortage(sum_from_dflow, sum_to_dflow)
+
+        # (self.demand["mf-riv2dflow1d_flux_negative"] - self.realised["dflow1d_flux2mf-riv_negative"]) * 60*60*24

imod_coupler/drivers/dfm_metamod/mf6_wrapper.py

@@ -175,15 +175,17 @@ def set_well_flux(
 
         if correction_flux is None or len(correction_flux) != len(flux):
             raise ValueError(f"Expected size of correction_flux is {len(flux)}")
-        for i in range(len(flux)):
-            flux[i, 0] = correction_flux[i]
+        flux[:, 0] = correction_flux[:]
+        # for i in range(len(flux)):
+        #     flux[i, 0] = correction_flux[i]
 
         self.set_value(bound_adress, flux)
 
-    def get_river_flux_estimate(
+    def get_river_drain_flux_estimate(
         self,
         mf6_flowmodel_key: str,
         mf6_river_pkg_key: str,
+        isdrain: Optional[bool] = False,
     ) -> NDArray[np.float_]:
         """
         returns the river1 fluxes consistent with current head, river stage and conductance.
@@ -221,8 +223,11 @@ def get_river_flux_estimate(
         nodelist = self.get_value_ptr(nodelist_adress)
 
         subset_head = head[nodelist - 1]
-        bot = bound[:, 2]
-        river_head = np.maximum(subset_head, bot)
+        if not isdrain:
+            bot = bound[:, 2]
+            river_head = np.maximum(subset_head, bot)
+        else:
+            river_head = subset_head
         q = NDArray[np.float_](len(nodelist))
         q[:] = bound[:, 1] * (bound[:, 0] - river_head)
 

tests/data/t_model_f/MODFLOW6/GWF_1/MODELINPUT/T-MODEL-F.NPF6

@@ -3,7 +3,6 @@
 #General Options
 
 BEGIN OPTIONS
- ALTERNATIVE_CELL_AVERAGING AMT-HMK
 END OPTIONS
 
 #Geology Options

tests/data/t_model_f/MODFLOW6/GWF_1/MODELINPUT/T-MODEL-F_SYS1.DRN6

@@ -13,6 +13,6 @@ MAXBOUND 2
 END DIMENSIONS
 
 BEGIN PERIOD 1
-    1     3     6    8.500000        100.0000         1
-    1     3     7    8.500000        100.0000         1
+    1     3     6    8.500000        100.0000
+    1     3     7    8.500000        100.0000
 END PERIOD