From a099533859b646af98d57e1142bc9ca23530f1fe Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Tue, 26 Mar 2024 15:39:11 +0100 Subject: [PATCH 1/7] Add newton classes --- imod_coupler/drivers/metamod/newton.py | 199 +++++++++++++++++++++++++ 1 file changed, 199 insertions(+) create mode 100644 imod_coupler/drivers/metamod/newton.py diff --git a/imod_coupler/drivers/metamod/newton.py b/imod_coupler/drivers/metamod/newton.py new file mode 100644 index 00000000..d93854c0 --- /dev/null +++ b/imod_coupler/drivers/metamod/newton.py @@ -0,0 +1,199 @@ +from typing import Any + +import numpy as np +from numpy.typing import NDArray + + +class Array: + """ + This class handles the MODFLOW 6 pointers and transforms them + from internal id's to user-id's and 3D shape (lay, row, col). + + ptr: pointer with shape (n_active_nodes) + flat: arrays with shape (n_user_nodes) + + """ + + ptr: NDArray[Any] + flat: NDArray[Any] + + def __init__( + self, + nlay: int, + nrow: int, + ncol: int, + userid: NDArray[Any], + ptr: NDArray[Any], + ) -> None: + self.nlay = nlay + self.nrow = nrow + self.ncol = ncol + self._set_ids(userid) + self.flat = np.zeros((self.nlay * self.nrow * self.ncol), dtype=np.float64) + self.ptr = ptr + + def update(self) -> None: + """ + updates flat array with pointer values + """ + self.flat[self.userid] = self.ptr + + def set_ptr(self, new_values: NDArray[Any]) -> None: + """ + sets pointer(n_active_nodes) with new values(n_user_nodes) + """ + valid = np.nonzero(new_values) + self.ptr[:] = new_values[valid] + + def _set_ids(self, userid: NDArray[Any]) -> None: + self.userid = userid - 1 + self.modelid = np.full((self.nlay * self.nrow * self.ncol), np.nan) + self.modelid[self.userid] = np.arange(self.userid.size) + + @property + def broadcasted(self) -> NDArray[Any]: + return self.flat.reshape((self.nlay, self.nrow, self.ncol)) + + @property + def sum(self) -> NDArray[Any]: + return np.sum(self.broadcasted, axis=0) + + +class BoundaryConditionArray(Array): + """ + This class handles the MODFLOW 6 pointers for boundary condition packages. The class handles + the extra transformation from list-oriented arrays to the model domain and back. + """ + + nodelist: NDArray[Any] + + def __init__( + self, + nlay: int, + nrow: int, + ncol: int, + userid: NDArray[Any], + ptr: NDArray[Any], + nodelist: NDArray[Any], + ) -> None: + self.nodelist = nodelist + super().__init__(nlay, nrow, ncol, userid, ptr) + + def update(self) -> None: + """ + updates flat array with package values. The flat array is of shape user_nodes. + The packages values are mapped using: + + 1- package nodeslist: package list to modeldomain(active_nodes) + 2- userid: modeldomain(active_nodes) to modeldomain(user_nodes) + + """ + model_index = self.nodelist - 1 + # remove invalid indices in case package nbound < maxbound + valid_index = model_index >= 0 + model_index = model_index[valid_index] + # zero array in case current nbound < previous nbound + self.flat[:] = 0.0 + self.flat[self.userid[model_index]] = self.ptr[:, 0][valid_index] + + def set_ptr(self, new_values: NDArray[Any]) -> None: + """ + Sets package pointers. This includes: + + 1- The package bounds array with boundary condition values + 2- The package nodeslist for mapping to modeldomain + + """ + valid = np.nonzero(new_values) + self.ptr[:, 0] = new_values[valid] + self.nodelist[:] = new_values[valid] + + +class PhreaticArray: + row_indices: NDArray[Any] + col_indices: NDArray[Any] + + def __init__( + self, + nlay: int, + nrow: int, + ncol: int, + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_package: NDArray[Any], + ) -> None: + self.nlay = nlay + self.nrow = nrow + self.ncol = ncol + self.row_indices = np.repeat(np.arange(nrow), ncol) + self.col_indices = np.tile(np.arange(ncol), reps=nrow) + self.model_condition = Array( + nlay=nlay, nrow=nrow, ncol=nrow, userid=userid, ptr=ptr_package + ) + self.saturation = Array( + nlay=nlay, nrow=nrow, ncol=nrow, userid=userid, ptr=ptr_saturation + ) + + def update(self) -> None: + self.model_condition.update() + self.saturation.update() + + @property + def phreatic_layer_indices(self) -> Any: + return np.argmax(self.saturation.broadcasted > 0, axis=0).flatten() + + @property + def phreatic_values(self) -> Any: + return self.model_condition.broadcasted[ + self.phreatic_layer_indices, self.row_indices, self.col_indices + ] + + +class PhreaticPackage(PhreaticArray): + def __init__( + self, + nlay: int, + nrow: int, + ncol: int, + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_package: NDArray[Any], + ) -> None: + self.nlay = nlay + self.nrow = nrow + self.ncol = ncol + super().__init__(nlay, nrow, ncol, userid, ptr_saturation, ptr_package) + + def set_pointer(self, array: NDArray[Any]) -> None: + self.model_condition.set_ptr(array) + + +class PhreaticBCPackage(PhreaticArray): + """ + This class assigns boundary condition values to the phreatic layer + based on the node saturation. The phreatic layer is defined at the top node + with a saturation > 0.0, along dimension layer. + """ + + def __init__( + self, + nlay: int, + nrow: int, + ncol: int, + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_bc: NDArray[Any], + ptr_bc_nodelist: NDArray[Any], + ) -> None: + self.nlay = nlay + self.nrow = nrow + self.ncol = ncol + super().__init__(nlay, nrow, ncol, userid, ptr_saturation, ptr_bc) + self.model_condition = BoundaryConditionArray( + nlay=nlay, + nrow=nrow, + ncol=nrow, + userid=userid, + ptr=ptr_bc, + nodelist=ptr_bc_nodelist, + ) From 14d40a7cff9b1cd9409787647dec27739184f241 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Thu, 2 May 2024 12:02:48 +0200 Subject: [PATCH 2/7] Add Newton kernel wrapper + create MetaMod instance for Newton formulation --- imod_coupler/config.py | 1 + imod_coupler/drivers/driver.py | 7 +- imod_coupler/drivers/metamod/metamod.py | 230 ++++++++++-- imod_coupler/drivers/metamod/newton.py | 199 ---------- .../drivers/ribametamod/ribametamod.py | 12 +- .../kernelwrappers/mf6_newton_wrapper.py | 339 ++++++++++++++++++ imod_coupler/kernelwrappers/mf6_wrapper.py | 31 +- .../primod/driver_coupling/metamod.py | 6 +- pre-processing/primod/metamod.py | 5 + tests/fixtures/fixture_metaswap.py | 24 +- tests/fixtures/fixture_modflow.py | 109 +++++- tests/test_imod_coupler/test_metamod.py | 121 +++++++ tests/test_imod_coupler/test_metamod_cases.py | 16 + tests/test_ribametamod.py | 216 ----------- 14 files changed, 859 insertions(+), 457 deletions(-) delete mode 100644 imod_coupler/drivers/metamod/newton.py create mode 100644 imod_coupler/kernelwrappers/mf6_newton_wrapper.py delete mode 100644 tests/test_ribametamod.py diff --git a/imod_coupler/config.py b/imod_coupler/config.py index d127fba7..acc54ab5 100644 --- a/imod_coupler/config.py +++ b/imod_coupler/config.py @@ -24,3 +24,4 @@ class BaseConfig(BaseModel): timing: bool = False driver_type: DriverType driver: BaseModel + modflow_newton_formulation: bool = False diff --git a/imod_coupler/drivers/driver.py b/imod_coupler/drivers/driver.py index dbe6480a..8779fffa 100644 --- a/imod_coupler/drivers/driver.py +++ b/imod_coupler/drivers/driver.py @@ -64,7 +64,7 @@ def get_driver( config_dict: dict[str, Any], config_dir: Path, base_config: BaseConfig ) -> Driver: from imod_coupler.drivers.metamod.config import MetaModConfig - from imod_coupler.drivers.metamod.metamod import MetaMod + from imod_coupler.drivers.metamod.metamod import MetaMod, MetaModNewton from imod_coupler.drivers.ribametamod.config import RibaMetaModConfig from imod_coupler.drivers.ribametamod.ribametamod import RibaMetaMod from imod_coupler.drivers.ribamod.config import RibaModConfig @@ -72,7 +72,10 @@ def get_driver( if base_config.driver_type == "metamod": metamod_config = MetaModConfig(config_dir=config_dir, **config_dict["driver"]) - return MetaMod(base_config, metamod_config) + if base_config.modflow_newton_formulation: + return MetaModNewton(base_config, metamod_config) + else: + return MetaMod(base_config, metamod_config) elif base_config.driver_type == "ribamod": ribamod_config = RibaModConfig(config_dir=config_dir, **config_dict["driver"]) return RibaMod(base_config, ribamod_config) diff --git a/imod_coupler/drivers/metamod/metamod.py b/imod_coupler/drivers/metamod/metamod.py index 9e1b1ca3..7a1141fc 100644 --- a/imod_coupler/drivers/metamod/metamod.py +++ b/imod_coupler/drivers/metamod/metamod.py @@ -17,6 +17,11 @@ from imod_coupler.config import BaseConfig from imod_coupler.drivers.driver import Driver from imod_coupler.drivers.metamod.config import Coupling, MetaModConfig +from imod_coupler.kernelwrappers.mf6_newton_wrapper import ( + PhreaticHeads, + PhreaticRecharge, + PhreaticStorage, +) from imod_coupler.kernelwrappers.mf6_wrapper import Mf6Wrapper from imod_coupler.kernelwrappers.msw_wrapper import MswWrapper from imod_coupler.logging.exchange_collector import ExchangeCollector @@ -39,7 +44,7 @@ class MetaMod(Driver): mf6_head: NDArray[Any] # the hydraulic head array in the coupled model mf6_recharge: NDArray[np.float64] # the coupled recharge array from the RCH package - mf6_storage: NDArray[Any] # the specific storage array (ss) + mf6_ss: NDArray[Any] # the specific storage array (ss) mf6_has_sc1: bool # when true, specific storage in mf6 is given as a storage coefficient (sc1) mf6_area: NDArray[Any] # cell area (size:nodes) mf6_top: NDArray[Any] # top of cell (size:nodes) @@ -60,6 +65,9 @@ class MetaMod(Driver): # dict. with mask arrays for msw=>mod coupling mask_msw2mod: dict[str, NDArray[Any]] = {} + node_idx: NDArray[Any] # mf6 nodes coupled to msw + msw_idx: NDArray[Any] # msw nodes coupled to MF6 + def __init__(self, base_config: BaseConfig, metamod_config: MetaModConfig): """Constructs the `MetaMod` object""" self.base_config = base_config @@ -106,7 +114,7 @@ def couple(self) -> None: self.mf6_recharge = self.mf6.get_recharge( self.coupling.mf6_model, self.coupling.mf6_msw_recharge_pkg ) - self.mf6_storage = self.mf6.get_storage(self.coupling.mf6_model) + self.mf6_ss = self.mf6.get_ss(self.coupling.mf6_model) self.mf6_has_sc1 = self.mf6.has_sc1(self.coupling.mf6_model) self.mf6_area = self.mf6.get_area(self.coupling.mf6_model) self.mf6_top = self.mf6.get_top(self.coupling.mf6_model) @@ -136,17 +144,19 @@ def couple(self) -> None: table_node2svat: NDArray[np.int32] = np.loadtxt( self.coupling.mf6_msw_node_map, dtype=np.int32, ndmin=2 ) - node_idx = table_node2svat[:, 0] - 1 - msw_idx = [ - svat_lookup[table_node2svat[ii, 1], table_node2svat[ii, 2]] - for ii in range(len(table_node2svat)) - ] + self.node_idx = table_node2svat[:, 0] - 1 + self.msw_idx = np.array( + [ + svat_lookup[table_node2svat[ii, 1], table_node2svat[ii, 2]] + for ii in range(len(table_node2svat)) + ] + ) self.map_msw2mod["storage"], self.mask_msw2mod["storage"] = create_mapping( - msw_idx, - node_idx, + self.msw_idx, + self.node_idx, self.msw_storage.size, - self.mf6_storage.size, + self.mf6_ss.size, "sum", ) @@ -167,8 +177,8 @@ def couple(self) -> None: self.map_msw2mod["storage"] = conversion_matrix * self.map_msw2mod["storage"] self.map_mod2msw["head"], self.mask_mod2msw["head"] = create_mapping( - node_idx, - msw_idx, + self.node_idx, + self.msw_idx, self.mf6_head.size, self.msw_head.size, "avg", @@ -178,13 +188,15 @@ def couple(self) -> None: self.coupling.mf6_msw_recharge_map, dtype=np.int32, ndmin=2 ) rch_idx = table_rch2svat[:, 0] - 1 - msw_idx = [ - svat_lookup[table_rch2svat[ii, 1], table_rch2svat[ii, 2]] - for ii in range(len(table_rch2svat)) - ] + self.msw_idx = np.array( + [ + svat_lookup[table_rch2svat[ii, 1], table_rch2svat[ii, 2]] + for ii in range(len(table_rch2svat)) + ] + ) self.map_msw2mod["recharge"], self.mask_msw2mod["recharge"] = create_mapping( - msw_idx, + self.msw_idx, rch_idx, self.msw_volume.size, self.mf6_recharge.size, @@ -205,16 +217,18 @@ def couple(self) -> None: ndmin=2, ) well_idx = table_well2svat[:, 0] - 1 - msw_idx = [ - svat_lookup[table_well2svat[ii, 1], table_well2svat[ii, 2]] - for ii in range(len(table_well2svat)) - ] + self.msw_idx = np.array( + [ + svat_lookup[table_well2svat[ii, 1], table_well2svat[ii, 2]] + for ii in range(len(table_well2svat)) + ] + ) ( self.map_msw2mod["sprinkling"], self.mask_msw2mod["sprinkling"], ) = create_mapping( - msw_idx, + self.msw_idx, well_idx, self.msw_volume.size, self.mf6_sprinkling_wells.size, @@ -239,6 +253,9 @@ def update(self) -> None: if has_converged: logger.debug(f"MF6-MSW converged in {kiter} iterations") break + if not has_converged: + logger.debug("MF6-MSW did not converge") + raise ValueError("MF6-MSW did not converge") self.mf6.finalize_solve(1) self.mf6.finalize_time_step() @@ -257,12 +274,12 @@ def get_end_time(self) -> float: def exchange_msw2mod(self) -> None: """Exchange Metaswap to Modflow""" - self.mf6_storage[:] = ( - self.mask_msw2mod["storage"][:] * self.mf6_storage[:] + self.mf6_ss[:] = ( + self.mask_msw2mod["storage"][:] * self.mf6_ss[:] + self.map_msw2mod["storage"].dot(self.msw_storage)[:] ) self.exchange_logger.log_exchange( - "mf6_storage", self.mf6_storage, self.get_current_time() + "mf6_storage", self.mf6_ss, self.get_current_time() ) self.exchange_logger.log_exchange( "msw_storage", self.msw_storage, self.get_current_time() @@ -306,3 +323,166 @@ def report_timing_totals(self) -> None: total_msw = self.msw.report_timing_totals() total = total_mf6 + total_msw logger.info(f"Total elapsed time in numerical kernels: {total:0.4f} seconds") + + +class MetaModNewton(MetaMod): + """ + MetaModNewton: the coupling between MetaSWAP and MODFLOW 6, for the Newton formulation of MODFLOW 6 + """ + + # TODO: check for fixed_cell in coupled rch-package + # TODO: check if nodes with well packages are confined + # TODO: check if storage formulation + npf-formulation are both set to convertible + + def __init__(self, base_config: BaseConfig, metamod_config: MetaModConfig): + super().__init__(base_config, metamod_config) + + def couple(self) -> None: + super().couple() + # Mapping for heads, SY and SS based on the top layer subset of MODFLOW 6 model arrays. + # The PhreaticHeads + PhreaticStorage methods get and set values to the corresponding + # phreatic nodes of potential deeper layers. + first_layer_nodes = self.get_first_layer_nodes() + self.map_mod2msw["head"], self.mask_mod2msw["head"] = create_mapping( + self.node_idx, + self.msw_idx, + first_layer_nodes.size, + self.msw_head.size, + "avg", + ) + + # mapping for SS now to top layer subset + self.map_msw2mod["storage"], self.mask_msw2mod["storage"] = create_mapping( + self.msw_idx, + self.node_idx, + self.msw_storage.size, + first_layer_nodes.size, + "sum", + ) + conversion_matrix = self.sto_conversion_terms( + self.mf6_has_sc1, first_layer_nodes + ) + self.map_msw2mod["storage"] = conversion_matrix * self.map_msw2mod["storage"] + # Create extra mapping for SY, since the SS mapping could contains a different conversion term + self.map_msw2mod["storage_sy"], self.mask_msw2mod["storage_sy"] = ( + create_mapping( + self.msw_idx, + self.node_idx, + self.msw_storage.size, + first_layer_nodes.size, + "sum", + ) + ) + # For exchange to SY, only the multiplication by area needs to be undone. + conversion_matrix = self.sto_conversion_terms(False, first_layer_nodes) + self.map_msw2mod["storage_sy"] = ( + conversion_matrix * self.map_msw2mod["storage_sy"] + ) + # Add Newton related phreatic exchange classes + self.couple_phreatic(first_layer_nodes) + + def couple_phreatic(self, first_layer_nodes: NDArray[Any]) -> None: + userid = self.mf6_get_userid() + saturation = self.mf6.get_saturation(self.coupling.mf6_model) + self.sy = self.mf6.get_sy(self.coupling.mf6_model) + self.sy_top = self.sy[first_layer_nodes] + self.ss_top = self.mf6_ss[first_layer_nodes] + self.recharge_nodelist = self.mf6.get_recharge_nodes( + self.coupling.mf6_model, self.coupling.mf6_msw_recharge_pkg + ) + self.hds = PhreaticHeads( + self.mf6.get_dis_shape(self.coupling.mf6_model), + userid, + saturation, + self.mf6_head, + first_layer_nodes, + ) + self.sto = PhreaticStorage( + self.mf6.get_dis_shape(self.coupling.mf6_model), + userid, + saturation, + self.sy, + self.mf6_ss, + first_layer_nodes, + ) + self.rch = PhreaticRecharge( + self.mf6.get_dis_shape(self.coupling.mf6_model), + userid, + saturation, + self.mf6_recharge, + self.recharge_nodelist, + ) + + def exchange_msw2mod(self) -> None: + """Exchange Metaswap to Modflow""" + new_sy = ( + self.mask_msw2mod["storage_sy"][:] * self.sy_top[:] + + self.map_msw2mod["storage_sy"].dot(self.msw_storage)[:] + ) + new_ss = ( + self.mask_msw2mod["storage"][:] * self.ss_top[:] + + self.map_msw2mod["storage"].dot(self.msw_storage)[:] + ) + self.sto.reset() + self.sto.set(new_sy, new_ss) + self.exchange_logger.log_exchange("mf6_sy", new_sy, self.get_current_time()) + self.exchange_logger.log_exchange( + "msw_storage", self.msw_storage, self.get_current_time() + ) + + # Set recharge + new_recharge = ( + self.mask_msw2mod["recharge"][:] * self.mf6_recharge[:] + + self.map_msw2mod["recharge"].dot(self.msw_volume)[:] / self.delt + ) / self.mf6_area[self.recharge_nodelist - 1] + self.rch.set(new_recharge) + + if self.enable_sprinkling_groundwater: + self.mf6_sprinkling_wells[:] = ( + self.mask_msw2mod["sprinkling"][:] * self.mf6_sprinkling_wells[:] + + self.map_msw2mod["sprinkling"].dot(self.msw_volume)[:] / self.delt + ) + + def exchange_mod2msw(self) -> None: + """Exchange Modflow to Metaswap""" + mf6_head = self.hds.get() + self.msw_head[:] = ( + self.mask_mod2msw["head"][:] * self.msw_head[:] + + self.map_mod2msw["head"].dot(mf6_head)[:] + ) + + def get_first_layer_nodes(self) -> NDArray[Any]: + _, nrow, ncol = self.mf6.get_dis_shape(self.coupling.mf6_model) + userid = self.mf6_get_userid() + first_layer_ids = userid[userid <= (nrow * ncol)] + if self.node_idx.max() > first_layer_ids.max(): + raise ValueError( + "MetaSWAP could only be coupled to the first model layer of MODFLOW 6" + ) + return first_layer_ids + + def mf6_get_userid(self) -> NDArray[Any]: + nlay, nrow, ncol = self.mf6.get_dis_shape(self.coupling.mf6_model) + userid = self.mf6.get_nodeuser(self.coupling.mf6_model) + if userid.size == 1: + # no reduced domain, set userid to modelid + # TODO: find out if there is a flag that indicates that usernodes == modelnodes + userid = np.arange(nlay * nrow * ncol) + 1 + return userid + + def sto_conversion_terms( + self, mf6_has_sc1: bool, indices: NDArray[Any] + ) -> NDArray[Any]: + if mf6_has_sc1: + conversion_terms = 1.0 / self.mf6_area[indices] + else: + conversion_terms = 1.0 / ( + self.mf6_area[indices] * (self.mf6_top[indices] - self.mf6_bot[indices]) + ) + + conversion_matrix = dia_matrix( + (conversion_terms, [0]), + shape=(self.mf6_area[indices].size, self.mf6_area[indices].size), + dtype=self.mf6_area.dtype, + ) + return conversion_matrix diff --git a/imod_coupler/drivers/metamod/newton.py b/imod_coupler/drivers/metamod/newton.py deleted file mode 100644 index d93854c0..00000000 --- a/imod_coupler/drivers/metamod/newton.py +++ /dev/null @@ -1,199 +0,0 @@ -from typing import Any - -import numpy as np -from numpy.typing import NDArray - - -class Array: - """ - This class handles the MODFLOW 6 pointers and transforms them - from internal id's to user-id's and 3D shape (lay, row, col). - - ptr: pointer with shape (n_active_nodes) - flat: arrays with shape (n_user_nodes) - - """ - - ptr: NDArray[Any] - flat: NDArray[Any] - - def __init__( - self, - nlay: int, - nrow: int, - ncol: int, - userid: NDArray[Any], - ptr: NDArray[Any], - ) -> None: - self.nlay = nlay - self.nrow = nrow - self.ncol = ncol - self._set_ids(userid) - self.flat = np.zeros((self.nlay * self.nrow * self.ncol), dtype=np.float64) - self.ptr = ptr - - def update(self) -> None: - """ - updates flat array with pointer values - """ - self.flat[self.userid] = self.ptr - - def set_ptr(self, new_values: NDArray[Any]) -> None: - """ - sets pointer(n_active_nodes) with new values(n_user_nodes) - """ - valid = np.nonzero(new_values) - self.ptr[:] = new_values[valid] - - def _set_ids(self, userid: NDArray[Any]) -> None: - self.userid = userid - 1 - self.modelid = np.full((self.nlay * self.nrow * self.ncol), np.nan) - self.modelid[self.userid] = np.arange(self.userid.size) - - @property - def broadcasted(self) -> NDArray[Any]: - return self.flat.reshape((self.nlay, self.nrow, self.ncol)) - - @property - def sum(self) -> NDArray[Any]: - return np.sum(self.broadcasted, axis=0) - - -class BoundaryConditionArray(Array): - """ - This class handles the MODFLOW 6 pointers for boundary condition packages. The class handles - the extra transformation from list-oriented arrays to the model domain and back. - """ - - nodelist: NDArray[Any] - - def __init__( - self, - nlay: int, - nrow: int, - ncol: int, - userid: NDArray[Any], - ptr: NDArray[Any], - nodelist: NDArray[Any], - ) -> None: - self.nodelist = nodelist - super().__init__(nlay, nrow, ncol, userid, ptr) - - def update(self) -> None: - """ - updates flat array with package values. The flat array is of shape user_nodes. - The packages values are mapped using: - - 1- package nodeslist: package list to modeldomain(active_nodes) - 2- userid: modeldomain(active_nodes) to modeldomain(user_nodes) - - """ - model_index = self.nodelist - 1 - # remove invalid indices in case package nbound < maxbound - valid_index = model_index >= 0 - model_index = model_index[valid_index] - # zero array in case current nbound < previous nbound - self.flat[:] = 0.0 - self.flat[self.userid[model_index]] = self.ptr[:, 0][valid_index] - - def set_ptr(self, new_values: NDArray[Any]) -> None: - """ - Sets package pointers. This includes: - - 1- The package bounds array with boundary condition values - 2- The package nodeslist for mapping to modeldomain - - """ - valid = np.nonzero(new_values) - self.ptr[:, 0] = new_values[valid] - self.nodelist[:] = new_values[valid] - - -class PhreaticArray: - row_indices: NDArray[Any] - col_indices: NDArray[Any] - - def __init__( - self, - nlay: int, - nrow: int, - ncol: int, - userid: NDArray[Any], - ptr_saturation: NDArray[Any], - ptr_package: NDArray[Any], - ) -> None: - self.nlay = nlay - self.nrow = nrow - self.ncol = ncol - self.row_indices = np.repeat(np.arange(nrow), ncol) - self.col_indices = np.tile(np.arange(ncol), reps=nrow) - self.model_condition = Array( - nlay=nlay, nrow=nrow, ncol=nrow, userid=userid, ptr=ptr_package - ) - self.saturation = Array( - nlay=nlay, nrow=nrow, ncol=nrow, userid=userid, ptr=ptr_saturation - ) - - def update(self) -> None: - self.model_condition.update() - self.saturation.update() - - @property - def phreatic_layer_indices(self) -> Any: - return np.argmax(self.saturation.broadcasted > 0, axis=0).flatten() - - @property - def phreatic_values(self) -> Any: - return self.model_condition.broadcasted[ - self.phreatic_layer_indices, self.row_indices, self.col_indices - ] - - -class PhreaticPackage(PhreaticArray): - def __init__( - self, - nlay: int, - nrow: int, - ncol: int, - userid: NDArray[Any], - ptr_saturation: NDArray[Any], - ptr_package: NDArray[Any], - ) -> None: - self.nlay = nlay - self.nrow = nrow - self.ncol = ncol - super().__init__(nlay, nrow, ncol, userid, ptr_saturation, ptr_package) - - def set_pointer(self, array: NDArray[Any]) -> None: - self.model_condition.set_ptr(array) - - -class PhreaticBCPackage(PhreaticArray): - """ - This class assigns boundary condition values to the phreatic layer - based on the node saturation. The phreatic layer is defined at the top node - with a saturation > 0.0, along dimension layer. - """ - - def __init__( - self, - nlay: int, - nrow: int, - ncol: int, - userid: NDArray[Any], - ptr_saturation: NDArray[Any], - ptr_bc: NDArray[Any], - ptr_bc_nodelist: NDArray[Any], - ) -> None: - self.nlay = nlay - self.nrow = nrow - self.ncol = ncol - super().__init__(nlay, nrow, ncol, userid, ptr_saturation, ptr_bc) - self.model_condition = BoundaryConditionArray( - nlay=nlay, - nrow=nrow, - ncol=nrow, - userid=userid, - ptr=ptr_bc, - nodelist=ptr_bc_nodelist, - ) diff --git a/imod_coupler/drivers/ribametamod/ribametamod.py b/imod_coupler/drivers/ribametamod/ribametamod.py index 271e79c5..42a80bfa 100644 --- a/imod_coupler/drivers/ribametamod/ribametamod.py +++ b/imod_coupler/drivers/ribametamod/ribametamod.py @@ -51,7 +51,7 @@ class RibaMetaMod(Driver): mf6_head: NDArray[Any] # the hydraulic head array in the coupled model mf6_recharge: NDArray[Any] # the coupled recharge array from the RCH package mf6_recharge_nodes: NDArray[Any] # node selection of rch nodes - mf6_storage: NDArray[Any] # the specific storage array (ss) + mf6_ss: NDArray[Any] # the specific storage array (ss) mf6_has_sc1: bool # when true, specific storage in mf6 is given as a storage coefficient (sc1) mf6_area: NDArray[Any] # cell area (size:nodes) mf6_top: NDArray[Any] # top of cell (size:nodes) @@ -220,7 +220,7 @@ def couple_metaswap(self) -> dict[str, Any]: self.mf6_recharge_nodes = self.mf6.get_recharge_nodes( self.coupling.mf6_model, self.coupling.mf6_msw_recharge_pkg ) - self.mf6_storage = self.mf6.get_storage(self.coupling.mf6_model) + self.mf6_ss = self.mf6.get_ss(self.coupling.mf6_model) self.mf6_has_sc1 = self.mf6.has_sc1(self.coupling.mf6_model) self.mf6_area = self.mf6.get_area(self.coupling.mf6_model) self.mf6_top = self.mf6.get_top(self.coupling.mf6_model) @@ -235,7 +235,7 @@ def couple_metaswap(self) -> dict[str, Any]: arrays["msw_storage"] = self.msw_storage arrays["mf6_recharge"] = self.mf6_recharge arrays["mf6_head"] = self.mf6_head - arrays["mf6_storage"] = self.mf6_storage + arrays["mf6_storage"] = self.mf6_ss arrays["mf6_has_sc1"] = self.mf6_has_sc1 arrays["mf6_area"] = self.mf6_area arrays["mf6_top"] = self.mf6_top @@ -444,12 +444,12 @@ def exchange_stage_rib2mod(self) -> None: def exchange_msw2mod(self) -> None: """Exchange Metaswap to Modflow""" - self.mf6_storage[:] = ( - self.mapping.msw2mod["storage_mask"][:] * self.mf6_storage[:] + self.mf6_ss[:] = ( + self.mapping.msw2mod["storage_mask"][:] * self.mf6_ss[:] + self.mapping.msw2mod["storage"].dot(self.msw_storage)[:] ) self.exchange_logger.log_exchange( - "mf6_storage", self.mf6_storage, self.get_current_time() + "mf6_storage", self.mf6_ss, self.get_current_time() ) self.exchange_logger.log_exchange( "msw_storage", self.msw_storage, self.get_current_time() diff --git a/imod_coupler/kernelwrappers/mf6_newton_wrapper.py b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py new file mode 100644 index 00000000..d626a23a --- /dev/null +++ b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py @@ -0,0 +1,339 @@ +from typing import Any + +import numpy as np +from numpy.typing import NDArray + + +class ExpandArray: + """ + This class handles the MODFLOW 6 internal arrays and transforms them + from internal size to user size and 3D shape (layer, row, col). + + Relative from the users input, the internal arrays are reduces in two ways: + 1- Reduction due to inactive model nodes (idomain != 1) + 2- Subset of model nodes for boundary condition package arrays. + + Parameters + ---------- + shape: list[int, int, int]: list that contains the grid demensions layer, row and column + userid (NDArray[Any]): array with user id's for model nodes + ptr (NDArray[Any]): pointer array with variable values on model nodes + ptr_nodelist (NDArray[Any]): optional pointer array for bc-packages to map package nodes to modelnodes + + """ + + variable_model: NDArray[Any] + variable_user: NDArray[Any] + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr: NDArray[Any], + ptr_nodelist: NDArray[Any] | None = None, + ) -> None: + self.nlay, self.nrow, self.ncol = shape + self._set_userid(userid) + self._set_modelid() + self.variable_model = ptr + self.variable_user = np.full( + (self.nlay * self.nrow * self.ncol), fill_value=np.nan, dtype=np.float64 + ) + self.nodelist = ptr_nodelist + + def _set_modelid(self) -> None: + """ + Sets array of 0-based modelid in user size + """ + self.modelid = np.full( + (self.nlay * self.nrow * self.ncol), np.nan, dtype=np.int32 + ) + self.modelid[self.userid] = np.arange(self.userid.size) + + def _set_userid(self, userid: NDArray[Any]) -> None: + """ + Sets array of 0-based userid in model size + """ + self.userid = userid - 1 + + def _update(self) -> None: + """ + updates variable array is user size with current pointer values from model sized array + """ + if self.nodelist is not None: + # boundary condition package; userid based on package nodelist, + # which is relative to the userid's + modelid = self.nodelist - 1 + self.variable_user[self.userid[modelid]] = self.variable_model[:] + else: + self.variable_user[self.userid] = self.variable_model[:] + + def broadcast(self) -> NDArray[Any]: + """ + returns 3D variable array is user size + """ + self._update() + return self.variable_user.reshape((self.nlay, self.nrow, self.ncol)) + + def expand(self) -> NDArray[Any]: + """ + returns flat variable array in user size + """ + self._update() + return self.variable_user + + +class PhreaticModelArray: + """ + This class set and gets phreatic elements from the pointer arrays for model-packages. + The initial nodelist gives the 'coupled' selection of nodes of the first model layer. + Based on the node saturation, the coresponding underlying phreatic nodes are selected + to get or set values. + + initial selection of nodes + | * | * | * | | + + heads phreatic nodes + layer 1 | -| | | | | * | | | | + layer 2 | | - | | | | | * | | | + layer 3 | | | - | - | | | | * | | + + Parameters + ---------- + shape: list[int, int, int]: list that contains the grid demensions layer, row and column + userid (NDArray[Any]): array with user id's for model nodes + ptr_saturation (NDArray[Any]): pointer array with saturation on model nodes + ptr_variable (NDArray[Any]): pointer array with variable values on model nodes + nodes (NDArray[Any]): selection of nodes in first model layer, for which the + corresponding phreaticnodes are computed + """ + + variable: ExpandArray + saturation: ExpandArray + nodes: NDArray[Any] + row_user_indices: NDArray[Any] + col_user_indices: NDArray[Any] + initialised: bool + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_variable: NDArray[Any], + nodes: NDArray[Any], + ) -> None: + self.nlay, self.nrow, self.ncol = shape + self.variable = ExpandArray(shape, userid, ptr_variable) + self.saturation = ExpandArray(shape, userid, ptr_saturation) + self.initialised = False + self.nodes = nodes - 1 + + def _set_user_indices(self) -> None: + """ + Computes the row and column indices for the user shape 3D array. For boundary + condition packages, the nodes array is only filled after the first prepare_timestep call + + considering inheritance, this method is therfore not called from the constructor + """ + self.row_user_indices = np.repeat(np.arange(self.nrow), self.ncol)[self.nodes] + self.col_user_indices = np.tile(np.arange(self.ncol), reps=self.nrow)[ + self.nodes + ] + self.initialised = True + + def set_ptr(self, new_values: NDArray[Any]) -> None: + """ + Sets the new values at the phreatic nodes + + Args: + new_values (NDArray[Any]): new values at nodes selection + """ + self.variable.variable_model[self.phreatic_modelid] = new_values + + def get_ptr(self) -> Any: + """ + Gets the variable values at the phreatic nodes + + Returns: + NDArray[Any]: Array with variable values at phreatic nodes of initial nodes selection + """ + return self.variable.broadcast()[ + self.phreatic_layer_user_indices, + self.row_user_indices, + self.col_user_indices, + ].flatten() # type ignore + + @property + def phreatic_layer_user_indices(self) -> Any: + self._set_user_indices() + return np.argmax(self.saturation.broadcast() > 0, axis=0).flatten()[self.nodes] + + @property + def phreatic_modelid(self) -> Any: + return self.variable.modelid.reshape((self.nlay, self.nrow, self.ncol))[ + self.phreatic_layer_user_indices, + self.row_user_indices, + self.col_user_indices, + ].flatten() + + +class PhreaticBCArray(PhreaticModelArray): + """ + This class sets the pointer arrays for the boundary condition-packages. In adition to the + PhreaticModelArray-class it sets the package nodelist pointer to the corresponing Phreatic nodes. + + This class set and gets phreatic elements from the pointer arrays for model-packages. + The initial nodelist gives the 'coupled' selection of nodes of the first model layer. + Based on the node saturation, the coresponding underlying phreatic nodes are selected + to get or set values. + """ + + initial_nodes: NDArray[Any] + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_package_variable: NDArray[Any], + ptr_package_nodelist: NDArray[Any], + ) -> None: + super().__init__( + shape, userid, ptr_saturation, ptr_package_variable, ptr_package_nodelist + ) + self.variable = ExpandArray( + shape, + userid, + ptr_package_variable, + ptr_package_nodelist, + ) + self.nodes_ptr = ptr_package_nodelist # used for nodes after prepare_timestep + + def _set_user_indices(self) -> None: + # the nodelist for bc-packages is only filled after the first prepare-timestep call + # This method therefore is not called from the constructor + if not self.initialised: + self.nodes = self.nodes_ptr - 1 + userid = self.variable.userid + self.row_user_indices = np.repeat(np.arange(self.nrow), self.ncol)[ + userid[self.nodes] + ] + self.col_user_indices = np.tile(np.arange(self.ncol), reps=self.nrow)[ + userid[self.nodes] + ] + self.initial_nodes = np.copy(self.nodes) + self.initialised = True + + def set_ptr(self, new_values: NDArray[Any]) -> None: + # The package nodelist is set to the phreatic nodes + self.variable.variable_model[:] = new_values[:] + self.variable.nodelist[:] = (self.phreatic_modelid + 1)[:] # type ignore + + @property + def phreatic_layer_user_indices(self) -> Any: + self._set_user_indices() + # use initial_nodes since self.nodes is updated by set_ptr method + return np.argmax(self.saturation.broadcast() > 0, axis=0).flatten()[ + self.initial_nodes + ] + + +class PhreaticStorage: + """ + Wrapper class for STO-package of MODFLOW 6. Contains methods for: + + 1- Setting values to phreatic nodes of SY and zeros the corresponding SS values + 2- Resetting both SY and SS arrays to initial values + """ + + sy: PhreaticModelArray + ss: PhreaticModelArray + zeros: NDArray[Any] + initial_sy: NDArray[Any] + initial_ss: NDArray[Any] + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_storage_sy: NDArray[Any], + ptr_storage_ss: NDArray[Any], + active_top_layer_nodes: NDArray[Any], + ) -> None: + self.zeros = np.zeros(active_top_layer_nodes.size) + self.initial_sy = np.copy(ptr_storage_sy) + self.initial_ss = np.copy(ptr_storage_ss) + self.sy = PhreaticModelArray( + shape, + userid, + ptr_saturation, + ptr_storage_sy, + active_top_layer_nodes, + ) + self.ss = PhreaticModelArray( + shape, + userid, + ptr_saturation, + ptr_storage_ss, + active_top_layer_nodes, + ) + + def reset(self) -> None: + self.sy.variable.variable_model[:] = self.initial_sy[:] + self.ss.variable.variable_model[:] = self.initial_ss[:] + + def set(self, new_sy: NDArray[Any], new_ss: NDArray[Any]) -> None: + self.sy.set_ptr(new_sy) + self.ss.set_ptr(self.zeros) + + +class PhreaticRecharge: + """ + Wrapper class for RCH-package of MODFLOW 6. Contains methods for: + + - Setting recharge values and updating package-nodelist to current phreatic nodes + """ + + recharge: PhreaticBCArray + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_recharge: NDArray[Any], + ptr_recharge_nodelist: NDArray[Any], + ) -> None: + self.recharge = PhreaticBCArray( + shape, userid, ptr_saturation, ptr_recharge, ptr_recharge_nodelist + ) + + def set(self, new_recharge: NDArray[Any]) -> None: + self.recharge.set_ptr(new_recharge) + + +class PhreaticHeads: + """ + Wrapper class for heads of MODFLOW 6. Contains methods for: + + - Getting heads values for current phreatic nodes + """ + + heads: PhreaticModelArray + + def __init__( + self, + shape: list[int], + userid: NDArray[Any], + ptr_saturation: NDArray[Any], + ptr_heads: NDArray[Any], + active_top_layer_nodes: NDArray[Any], + ) -> None: + self.heads = PhreaticModelArray( + shape, userid, ptr_saturation, ptr_heads, active_top_layer_nodes + ) + + def get(self) -> NDArray[Any]: + return self.heads.get_ptr() diff --git a/imod_coupler/kernelwrappers/mf6_wrapper.py b/imod_coupler/kernelwrappers/mf6_wrapper.py index e8d8fb4b..3b194773 100644 --- a/imod_coupler/kernelwrappers/mf6_wrapper.py +++ b/imod_coupler/kernelwrappers/mf6_wrapper.py @@ -40,6 +40,10 @@ def get_drainage_packages( key: Mf6Drainage(self, mf6_flowmodel_key, key) for key in mf6_drainage_keys } + def has_newton(self, mf6_flowmodel_key: str) -> bool: + newton_tag = self.get_var_address("INEWTON", mf6_flowmodel_key) + return self.get_value_ptr(newton_tag)[0] == 1 + def get_well( self, mf6_flowmodel_key: str, @@ -68,11 +72,20 @@ def get_recharge_nodes( ) return self.get_value_ptr(mf6_recharge_nodes_tag) - def get_storage(self, mf6_flowmodel_key: str) -> NDArray[np.float64]: + def get_ss(self, mf6_flowmodel_key: str) -> NDArray[np.float64]: mf6_storage_tag = self.get_var_address("SS", mf6_flowmodel_key, "STO") mf6_storage = self.get_value_ptr(mf6_storage_tag) return mf6_storage + def get_sy(self, mf6_flowmodel_key: str) -> NDArray[np.float64]: + mf6_storage_tag = self.get_var_address("SY", mf6_flowmodel_key, "STO") + mf6_storage = self.get_value_ptr(mf6_storage_tag) + return mf6_storage + + def get_saturation(self, mf6_flowmodel_key: str) -> NDArray[np.float64]: + saturation_tag = self.get_var_address("SAT", mf6_flowmodel_key, "NPF") + return self.get_value_ptr(saturation_tag) + def has_sc1(self, mf6_flowmodel_key: str) -> bool: mf6_is_sc1_tag = self.get_var_address("ISTOR_COEF", mf6_flowmodel_key, "STO") mf6_has_sc1 = bool(self.get_value_ptr(mf6_is_sc1_tag)[0] != 0) @@ -93,6 +106,20 @@ def get_bot(self, mf6_flowmodel_key: str) -> NDArray[np.float64]: mf6_bot = self.get_value_ptr(mf6_bot_tag) return mf6_bot + def get_dis_shape(self, mf6_flowmodel_key: str) -> tuple: + nlay_tag = self.get_var_address("NLAY", mf6_flowmodel_key, "DIS") + nrow_tag = self.get_var_address("NROW", mf6_flowmodel_key, "DIS") + ncol_tag = self.get_var_address("NCOL", mf6_flowmodel_key, "DIS") + return ( + self.get_value_ptr(nlay_tag)[0], + self.get_value_ptr(nrow_tag)[0], + self.get_value_ptr(ncol_tag)[0], + ) + + def get_nodeuser(self, mf6_flowmodel_key: str) -> NDArray[np.int32]: + nodeuser_tag = self.get_var_address("NODEUSER", mf6_flowmodel_key, "DIS") + return self.get_value_ptr(nodeuser_tag) + def max_iter(self) -> Any: mf6_max_iter_tag = self.get_var_address("MXITER", "SLN_1") mf6_max_iter = self.get_value_ptr(mf6_max_iter_tag)[0] @@ -126,7 +153,7 @@ def get_drainage_elevation( Returns ------- - NDArray[np.float64]: + NDArray[np.float64]: Drainage elevation in modflow """ bound_address = self.get_var_address( diff --git a/pre-processing/primod/driver_coupling/metamod.py b/pre-processing/primod/driver_coupling/metamod.py index 46dfe158..a49109e3 100644 --- a/pre-processing/primod/driver_coupling/metamod.py +++ b/pre-processing/primod/driver_coupling/metamod.py @@ -1,7 +1,7 @@ from pathlib import Path from typing import Any -from imod.mf6 import GroundwaterFlowModel +from imod.mf6 import GroundwaterFlowModel, Modflow6Simulation from imod.msw import GridData, MetaSwapModel, Sprinkling from primod.driver_coupling.driver_coupling_base import DriverCoupling @@ -27,6 +27,9 @@ class MetaModDriverCoupling(DriverCoupling): mf6_recharge_package: str mf6_wel_package: str | None = None + def has_newton_formulation(self, mf6_simulation: Modflow6Simulation) -> bool: + return bool(mf6_simulation[self.mf6_model]._options["newton"]) + def _check_sprinkling( self, msw_model: MetaSwapModel, gwf_model: GroundwaterFlowModel ) -> bool: @@ -96,7 +99,6 @@ def write_exchanges(self, directory: Path, coupled_model: Any) -> dict[str, Any] coupling_dict: dict[str, Any] = {} coupling_dict["mf6_model"] = self.mf6_model - coupling_dict["mf6_msw_node_map"] = grid_mapping.write(directory) coupling_dict["mf6_msw_recharge_pkg"] = self.mf6_recharge_package coupling_dict["mf6_msw_recharge_map"] = rch_mapping.write(directory) diff --git a/pre-processing/primod/metamod.py b/pre-processing/primod/metamod.py index 07783e94..60bcbdd4 100644 --- a/pre-processing/primod/metamod.py +++ b/pre-processing/primod/metamod.py @@ -34,6 +34,9 @@ def __init__( self.msw_model = msw_model self.mf6_simulation = mf6_simulation self.coupling_list = coupling_list + self.newton_formulation = self.coupling_list[0].has_newton_formulation( + self.mf6_simulation + ) def write( self, @@ -154,6 +157,8 @@ def write_toml( "coupling": [coupling_dict], }, } + if self.newton_formulation: + coupler_toml["modflow_newton_formulation"] = True with open(toml_path, "wb") as f: tomli_w.dump(coupler_toml, f) diff --git a/tests/fixtures/fixture_metaswap.py b/tests/fixtures/fixture_metaswap.py index 908b88e8..37153add 100644 --- a/tests/fixtures/fixture_metaswap.py +++ b/tests/fixtures/fixture_metaswap.py @@ -67,7 +67,7 @@ def metaswap_model( area, xr.full_like(area, 1, dtype=int), xr.full_like(area, 1.0, dtype=float), - xr.full_like(active, 1.0, dtype=float), + xr.full_like(active, 0.0, dtype=float), # to be consistent with mf6 model xr.full_like(active, 1, dtype=int), active, ) @@ -238,3 +238,25 @@ def prepared_msw_model_inactive( ) return msw_model + + +@pytest_cases.fixture(scope="function") +def prepared_msw_model_newton( + partly_inactive_idomain: xr.DataArray, + metaswap_lookup_table: Path, +) -> msw.MetaSwapModel: + msw_model = make_msw_model(partly_inactive_idomain) + # increase precipitation, zero evaporation + msw_model["meteo_grid"].dataset["precipitation"] = ( + msw_model["meteo_grid"].dataset["precipitation"] * 6.0 + ) + msw_model["meteo_grid"].dataset["evapotranspiration"] = ( + msw_model["meteo_grid"].dataset["evapotranspiration"] * 0.0 + ) + + # Override unsat_svat_path with path from environment + msw_model.simulation_settings["unsa_svat_path"] = ( + msw_model._render_unsaturated_database_path(metaswap_lookup_table) + ) + + return msw_model diff --git a/tests/fixtures/fixture_modflow.py b/tests/fixtures/fixture_modflow.py index 82369a82..d2e65b78 100644 --- a/tests/fixtures/fixture_modflow.py +++ b/tests/fixtures/fixture_modflow.py @@ -7,14 +7,16 @@ from .common import create_wells, get_times, grid_sizes -def make_mf6_model(idomain: xr.DataArray) -> mf6.GroundwaterFlowModel: +def make_mf6_model( + idomain: xr.DataArray, newton: bool = False +) -> mf6.GroundwaterFlowModel: _, _, layer, _, _, dz = grid_sizes() nlay = len(layer) top = 0.0 bottom = top - xr.DataArray(np.cumsum(dz), coords={"layer": layer}, dims="layer") - gwf_model = mf6.GroundwaterFlowModel() + gwf_model = mf6.GroundwaterFlowModel(newton=newton) gwf_model["dis"] = mf6.StructuredDiscretization( idomain=idomain, top=top, bottom=bottom ) @@ -72,10 +74,79 @@ def make_recharge_pkg(idomain: xr.DataArray) -> mf6.Recharge: recharge[:, 0] = np.nan recharge = recharge.where(idomain_l1) - return mf6.Recharge(recharge) + return mf6.Recharge(recharge, save_flows=True, fixed_cell=1) + + +def make_coupled_mf6_model_newton(idomain: xr.DataArray) -> mf6.Modflow6Simulation: + nlay, nrow, ncol = idomain.shape + gwf_model = make_mf6_model(idomain, newton=True) + # Add coupling package + gwf_model["rch_msw"] = make_recharge_pkg(idomain) + # Add GHB-package + times = get_times() + head = xr.full_like(idomain.astype(np.float64), np.nan) + head[-1, :, :] = -2.5 + head = head.expand_dims(time=times) + head = head.where(head.time > head.time[100]) + conductance = (xr.ones_like(head) * 1000).where(head.notnull()) + gwf_model["ghb"] = mf6.GeneralHeadBoundary( + head=head.where(idomain == 1), + conductance=conductance.where(idomain == 1), + print_input=True, + print_flows=True, + save_flows=True, + ) + # update dis-package + gwf_model.pop("dis") + top = 0.0 + dz = np.array([1.0] * nlay) + bottom = top - xr.DataArray( + np.cumsum(dz), coords={"layer": idomain.layer}, dims="layer" + ) + gwf_model["dis"] = mf6.StructuredDiscretization( + idomain=idomain, top=top, bottom=bottom + ) + # update npf-package + gwf_model.pop("npf") + k_values = np.ones(nlay) + k = xr.DataArray(k_values, {"layer": idomain.layer}, ("layer",)) + k33 = k + icelltype = xr.full_like(bottom, 1, dtype=int) + gwf_model["npf"] = mf6.NodePropertyFlow( + icelltype=icelltype, k=k, k33=k33, save_saturation=True, save_flows=True + ) + # update ic-package + gwf_model.pop("ic") + gwf_model["ic"] = mf6.InitialConditions(start=-2.5) + # sto-package + gwf_model["sto"].dataset["convertible"] = 1 + + # add coupling packages + gwf_model["rch_msw"] = make_recharge_pkg(idomain) + # gwf_model["wells_msw"] = create_wells(nrow, ncol, idomain) + + simulation = make_mf6_simulation(gwf_model) + simulation.pop("solver") + simulation["solver"] = mf6.Solution( + ["GWF_1"], + complexity="complex", + outer_dvclose=1e-5, + outer_maximum=500, + backtracking_number=0, + inner_maximum=100, + inner_dvclose=1e-6, + inner_rclose=1e-6, + rclose_option="strict", + linear_acceleration="bicgstab", + number_orthogonalizations=10, + relaxation_factor=0.0, + ) + return simulation -def make_coupled_mf6_model(idomain: xr.DataArray) -> mf6.Modflow6Simulation: +def make_coupled_mf6_model( + idomain: xr.DataArray, newton: bool = False +) -> mf6.Modflow6Simulation: _, nrow, ncol = idomain.shape gwf_model = make_mf6_model(idomain) times = get_times() @@ -133,6 +204,21 @@ def make_idomain() -> xr.DataArray: ) +def make_idomain_partly_inactive() -> xr.DataArray: + x, y, layer, dx, dy, _ = grid_sizes() + nlay = len(layer) + nrow = len(y) + ncol = len(x) + array = np.ones((nlay, nrow, ncol), dtype=np.int32) + array[:, :, 0] = 0 + + return xr.DataArray( + data=array, + dims=("layer", "y", "x"), + coords={"layer": layer, "y": y, "x": x, "dx": dx, "dy": dy}, + ) + + @pytest_cases.fixture(scope="function") def active_idomain() -> xr.DataArray: """Return all active idomain""" @@ -141,6 +227,14 @@ def active_idomain() -> xr.DataArray: return idomain +@pytest_cases.fixture(scope="function") +def partly_inactive_idomain() -> xr.DataArray: + """Return partly active idomain""" + idomain = make_idomain_partly_inactive() + + return idomain + + @pytest_cases.fixture(scope="function") def inactive_idomain() -> xr.DataArray: """Return idomain with an inactive cell""" @@ -161,6 +255,13 @@ def coupled_mf6_model(active_idomain: xr.DataArray) -> mf6.Modflow6Simulation: return make_coupled_mf6_model(active_idomain) +@pytest_cases.fixture(scope="function") +def coupled_mf6_model_newton( + partly_inactive_idomain: xr.DataArray, +) -> mf6.Modflow6Simulation: + return make_coupled_mf6_model_newton(partly_inactive_idomain) + + @pytest_cases.fixture(scope="function") def mf6_bucket_model(active_idomain: xr.DataArray) -> mf6.Modflow6Simulation: # The bottom of the ribasim trivial model is located at 0.0 m: the surface diff --git a/tests/test_imod_coupler/test_metamod.py b/tests/test_imod_coupler/test_metamod.py index ac7d343b..b097a2c7 100644 --- a/tests/test_imod_coupler/test_metamod.py +++ b/tests/test_imod_coupler/test_metamod.py @@ -4,6 +4,7 @@ from collections.abc import Callable from pathlib import Path +import numpy as np import pytest import tomli import tomli_w @@ -14,6 +15,12 @@ from pytest_cases import parametrize_with_cases from test_utilities import numeric_csvfiles_equal +from imod_coupler.kernelwrappers.mf6_newton_wrapper import ( + PhreaticHeads, + PhreaticRecharge, + PhreaticStorage, +) + def mf6_output_files(path: Path) -> tuple[Path, Path, Path, Path]: """return paths to Modflow 6 output files""" @@ -146,6 +153,7 @@ def test_metamod_develop( modflow6_dll=modflow_dll_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 / metamod_model._toml_name) @@ -411,3 +419,116 @@ def add_logging_request_to_toml_file(toml_dir: Path, toml_filename: str) -> None ) # on print ,"\\\\" gets rendered as "\\" with open(toml_dir / "output_config.toml", "w") as f: f.write(output_config_content.format(workdir=path_quadruple_backslash)) + + +def test_newton_clases() -> None: + nlay = 3 + nrow = 4 + ncol = 4 + idomain = np.ones((nlay, nrow, ncol)) + idomain[1, :, :] = -1 # layer two is inactive + userid = (np.arange(nlay * nrow * ncol) + 1).reshape((nlay, nrow, ncol))[ + idomain == 1 + ] + recharge = np.full((3 * 4), fill_value=0.001) + recharge_nodelist = np.array( + [1, 2, 3, 5, 6, 7, 9, 10, 11, 13, 14, 15] + ) # layer 1, first 3 columns + + saturation = np.ones_like(idomain) + saturation[2, :, 2] = 0.5 # third column, phreatisch layer = 3 + saturation[0:2, :, 2] = 0.0 + saturation[0, :, 1] = 0.6 # second column, phreatisch layer = 1 + saturation[0, :, 0] = 1.0 # first column, phreatisch layer = 1 + saturation = saturation[idomain == 1] + + phreatic_nodelist = [ + 1, + 2, + 35 - 16, + 5, + 6, + 39 - 16, + 9, + 10, + 43 - 16, + 13, + 14, + 47 - 16, + ] # -16 for userid -> modelid since layer 2 is inactive + + rch = PhreaticRecharge( + [nlay, nrow, ncol], userid, saturation, recharge, recharge_nodelist + ) + recharge_org = np.copy(recharge) + rch.set(recharge * 2) + + assert (recharge == recharge_org * 2).all() + assert (recharge_nodelist == phreatic_nodelist).all() + + # storage + sy = np.full( + (nlay - 1, nrow, ncol), fill_value=0.15 + ).flatten() # -1 to remove inactive second layer from internal array + sy_org = np.copy(sy) + ss = np.full( + (nlay - 1, nrow, ncol), fill_value=0.1e-6 + ).flatten() # -1 to remove inactive second layer from internal array + ss_org = np.copy(ss) + coupled_nodes = np.array([1, 2, 3, 5, 6, 7, 9, 10, 11, 13, 14, 15]) + + # set storage pointers based on saturation + sto = PhreaticStorage( + [nlay, nrow, ncol], + userid, + saturation, + sy, + ss, + coupled_nodes, + ) + + new_sy = np.full((coupled_nodes.size), fill_value=0.20) + # coupled first layer, first three columns + # should only updates the phreatic nodes underlying the coupled nodes + sto.set(new_sy, new_sy) + + phreatic_nodes = np.array( + [ + 1, + 2, + 35 - 16, + 5, + 6, + 39 - 16, + 9, + 10, + 43 - 16, + 13, + 14, + 47 - 16, + ] + ) # -16 for userid -> modelid since layer 2 is inactive + + nodes = np.arange((nlay - 1) * nrow * ncol) + 1 + mask = np.ones_like(nodes) + mask[phreatic_nodes - 1] = 0 + non_phreatic_nodes = nodes[mask.astype(dtype=bool)] + + # sets phreatic values for SY + assert (sy[phreatic_nodes - 1] == 0.2).all() + assert (sy[non_phreatic_nodes - 1] == 0.15).all() + # zeros SS for phreatic nodes + assert (ss[phreatic_nodes - 1] == 0.0).all() + assert (ss[non_phreatic_nodes - 1] == 0.1e-6).all() + + # resets arrays to initial values + sto.reset() + assert (sy == sy_org).all() + assert (ss == ss_org).all() + + # head + heads = np.arange((nlay - 1) * nrow * ncol) + + hds = PhreaticHeads([nlay, nrow, ncol], userid, saturation, heads, coupled_nodes) + phreatic_heads = hds.get() + assert (phreatic_heads == heads[phreatic_nodes - 1]).all() diff --git a/tests/test_imod_coupler/test_metamod_cases.py b/tests/test_imod_coupler/test_metamod_cases.py index 886a0338..dcfb4fae 100644 --- a/tests/test_imod_coupler/test_metamod_cases.py +++ b/tests/test_imod_coupler/test_metamod_cases.py @@ -190,3 +190,19 @@ def cases_metamod_sprinkling( ) return metamod_ss, metamod_sc + + +def case_newton( + coupled_mf6_model_newton: Modflow6Simulation, + prepared_msw_model_newton: MetaSwapModel, +) -> MetaMod: + prepared_msw_model_newton.pop("sprinkling") + + driver_coupling = MetaModDriverCoupling( + mf6_model="GWF_1", mf6_recharge_package="rch_msw" + ) + return MetaMod( + prepared_msw_model_newton, + coupled_mf6_model_newton, + coupling_list=[driver_coupling], + ) diff --git a/tests/test_ribametamod.py b/tests/test_ribametamod.py deleted file mode 100644 index 7accea23..00000000 --- a/tests/test_ribametamod.py +++ /dev/null @@ -1,216 +0,0 @@ -from collections.abc import Callable -from pathlib import Path - -import numpy as np -import pandas as pd -import pytest -from imod.msw import MetaSwapModel -from primod.ribametamod import RibaMetaMod -from pytest_cases import parametrize_with_cases - -from imod_coupler.drivers.ribametamod.exchange import ExchangeBalance - -msw_outputlabel_ponding: str = " Pssw(m3)" -msw_outputlabel_swsprinkling: str = " ts2dfmput(m3)" - - -@pytest.mark.xfail(reason="MetaSWAP issues") -@pytest.mark.xdist_group(name="ribasim") -@parametrize_with_cases("ribametamod_model", glob="backwater_model") -def test_ribametamod_backwater( - tmp_path_dev: Path, - ribametamod_model: RibaMetaMod | MetaSwapModel, - 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], - ribametamod_backwater_tot_svat_ref: Path, -) -> None: - """ - Test if the backwater model works as expected - """ - 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, - ) - run_coupler_function(tmp_path_dev / ribametamod_model._toml_name) - tot_svat_reference = pd.read_csv(ribametamod_backwater_tot_svat_ref) - tot_svat_test = pd.read_csv("metaswap/msw/csv/tot_svat_per.csv") - assert tot_svat_test[msw_outputlabel_swsprinkling].equals( - tot_svat_reference[msw_outputlabel_swsprinkling] - ) - assert tot_svat_test[msw_outputlabel_ponding].equals( - tot_svat_reference[msw_outputlabel_ponding] - ) - - -@pytest.mark.xfail(reason="MetaSWAP issues") -@pytest.mark.xdist_group(name="ribasim") -@parametrize_with_cases("ribametamod_model", glob="bucket_model") -def test_ribametamod_bucket( - tmp_path_dev: Path, - ribametamod_model: RibaMetaMod | MetaSwapModel, - 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], - ribametamod_bucket_tot_svat_ref: Path, -) -> None: - """ - Test if the bucket model works as expected - """ - 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, - ) - run_coupler_function(tmp_path_dev / ribametamod_model._toml_name) - tot_svat_reference = pd.read_csv(ribametamod_bucket_tot_svat_ref) - tot_svat_test = pd.read_csv("metaswap/msw/csv/tot_svat_per.csv") - assert tot_svat_test[msw_outputlabel_swsprinkling].equals( - tot_svat_reference[msw_outputlabel_swsprinkling] - ) - assert tot_svat_test[msw_outputlabel_ponding].equals( - tot_svat_reference[msw_outputlabel_ponding] - ) - - -@pytest.mark.xfail(reason="MetaSWAP issues") -@pytest.mark.xdist_group(name="ribasim") -@parametrize_with_cases("ribametamod_model", glob="two_basin_model") -def test_ribametamod_two_basin( - tmp_path_dev: Path, - ribametamod_model: RibaMetaMod, - metaswap_dll_devel: Path, - metaswap_dll_dep_dir_devel: Path, - modflow_dll_devel: Path, - ribasim_dll_devel: Path, - ribasim_dll_dep_dir_devel: Path, - run_coupler_function: Callable[[Path], None], - ribametamod_two_basin_tot_svat_ref: Path, -) -> None: - """ - Test if the two-basin model model works as expected - """ - 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, - ) - - run_coupler_function(tmp_path_dev / ribametamod_model._toml_name) - tot_svat_reference = pd.read_csv(ribametamod_two_basin_tot_svat_ref) - tot_svat_test = pd.read_csv("metaswap/msw/csv/tot_svat_per.csv") - assert tot_svat_test[msw_outputlabel_swsprinkling].equals( - tot_svat_reference[msw_outputlabel_swsprinkling] - ) - assert tot_svat_test[msw_outputlabel_ponding].equals( - tot_svat_reference[msw_outputlabel_ponding] - ) - - -def test_exchange_balance() -> None: - shape = 4 - labels = ["flux-1", "flux-2"] - exchange = ExchangeBalance(shape=shape, labels=labels) - - # exchange demands to class - array_negative = np.zeros(shape=shape, dtype=np.float64) - array_positive = np.zeros(shape=shape, dtype=np.float64) - - # seperate negative contributions for n:1 exchange - array_negative[0] = -10 - array_negative[1] = -10 - array_positive[0] = 0.0 - array_positive[1] = 5.0 - - demand_array = array_negative + array_positive - exchange.demands["flux-1"] = demand_array - exchange.demands["flux-2"] = demand_array * 0.5 - exchange.demands_negative["flux-1"] = array_negative - exchange.demands_negative["flux-2"] = array_negative * 0.5 - - # check summed demand - assert np.all(exchange.demand == demand_array + (demand_array * 0.5)) - # check summed negative demand - assert np.all(exchange.demand_negative == array_negative + (array_negative * 0.5)) - - # evaluate realised method - realised = np.zeros(shape=shape, dtype=np.float64) - realised[0] = -5.0 - realised[1] = -5.0 - # compute - exchange.compute_realised(realised) - # compare: realised_factor = 1 - (-shortage - sum_negative_demands) - realised_factor = np.zeros(shape=shape, dtype=np.float64) - realised_factor[0] = 1 - (-10 / -15) - realised_factor[1] = 1 - (-2.5 / -15) - - expected_flux1 = np.zeros(shape=shape, dtype=np.float64) - expected_flux2 = np.zeros(shape=shape, dtype=np.float64) - expected_flux1[0] = realised_factor[0] * array_negative[0] - expected_flux2[0] = realised_factor[0] * array_negative[0] * 0.5 - expected_flux1[1] = realised_factor[1] * array_negative[1] - expected_flux2[1] = realised_factor[1] * array_negative[1] * 0.5 - assert np.all(expected_flux1 == exchange.realised_negative["flux-1"]) - assert np.all(expected_flux2 == exchange.realised_negative["flux-2"]) - - compute_realised = np.zeros(shape=shape, dtype=np.float64) - compute_realised[0] = ( - exchange.realised_negative["flux-1"][0] - + exchange.realised_negative["flux-2"][0] - + array_positive[0] - + (array_positive[0] * 0.5) - ) - compute_realised[1] = ( - exchange.realised_negative["flux-1"][1] - + exchange.realised_negative["flux-2"][1] - + array_positive[1] - + (array_positive[1] * 0.5) - ) - assert np.all(np.isclose(realised, compute_realised)) - - # check if reset zeros arrays - exchange.reset() - assert np.all(exchange.demand == np.zeros(shape=shape, dtype=np.float64)) - assert np.all(exchange.demand_negative == np.zeros(shape=shape, dtype=np.float64)) - - # check if errors are thrown - # shortage larger than negative demands - shape = 1 - labels = ["flux-1"] - exchange = ExchangeBalance(shape=shape, labels=labels) - exchange.demands["flux-1"] = np.ones(shape=shape, dtype=np.float64) * -4 - exchange.demands_negative["flux-1"] = np.ones(shape=shape, dtype=np.float64) * -4 - realised = np.ones(shape=shape, dtype=np.float64) * 2 - with pytest.raises( - ValueError, - match="Invalid realised values: found shortage larger than negative demand contributions", - ): - exchange.compute_realised(realised) - - # shortage for positive demand - shape = 1 - labels = ["flux-1"] - exchange = ExchangeBalance(shape=shape, labels=labels) - exchange.demands["flux-1"] = np.ones(shape=shape, dtype=np.float64) * 10 - exchange.demands_negative["flux-1"] = np.ones(shape=shape, dtype=np.float64) * -4 - realised = np.ones(shape=shape, dtype=np.float64) * 8 - with pytest.raises( - ValueError, match="Invalid realised values: found shortage for positive demand" - ): - exchange.compute_realised(realised) From 1aa0336a1d0cfe6d124742b5f03f8d61b882689d Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Thu, 2 May 2024 17:17:13 +0200 Subject: [PATCH 3/7] Add checks to primod --- imod_coupler/drivers/metamod/metamod.py | 4 -- .../primod/driver_coupling/metamod.py | 55 ++++++++++++++++++- tests/fixtures/fixture_modflow.py | 4 +- 3 files changed, 55 insertions(+), 8 deletions(-) diff --git a/imod_coupler/drivers/metamod/metamod.py b/imod_coupler/drivers/metamod/metamod.py index 7a1141fc..56bb2a3a 100644 --- a/imod_coupler/drivers/metamod/metamod.py +++ b/imod_coupler/drivers/metamod/metamod.py @@ -330,10 +330,6 @@ class MetaModNewton(MetaMod): MetaModNewton: the coupling between MetaSWAP and MODFLOW 6, for the Newton formulation of MODFLOW 6 """ - # TODO: check for fixed_cell in coupled rch-package - # TODO: check if nodes with well packages are confined - # TODO: check if storage formulation + npf-formulation are both set to convertible - def __init__(self, base_config: BaseConfig, metamod_config: MetaModConfig): super().__init__(base_config, metamod_config) diff --git a/pre-processing/primod/driver_coupling/metamod.py b/pre-processing/primod/driver_coupling/metamod.py index a49109e3..dea7c492 100644 --- a/pre-processing/primod/driver_coupling/metamod.py +++ b/pre-processing/primod/driver_coupling/metamod.py @@ -1,6 +1,8 @@ from pathlib import Path from typing import Any +import xarray as xr +from imod import mf6 from imod.mf6 import GroundwaterFlowModel, Modflow6Simulation from imod.msw import GridData, MetaSwapModel, Sprinkling @@ -28,7 +30,51 @@ class MetaModDriverCoupling(DriverCoupling): mf6_wel_package: str | None = None def has_newton_formulation(self, mf6_simulation: Modflow6Simulation) -> bool: - return bool(mf6_simulation[self.mf6_model]._options["newton"]) + has_newton = bool(mf6_simulation[self.mf6_model]._options["newton"]) + if has_newton: + self._check_newton_simulation_settings(mf6_simulation[self.mf6_model]) + return has_newton + + def _check_newton_simulation_settings( + self, gwf_model: GroundwaterFlowModel + ) -> None: + # check if both npf-package and sto-package are convertible since npf-sat is used in coupling + # and SY is used as exchange variable + idomain = get_idomain(gwf_model) + mask = xr.ones_like(idomain) + for label in gwf_model: + package = gwf_model[label] + if isinstance(package, mf6.NodePropertyFlow): + # broadcast (layered) constants + npf_celtype = (mask * gwf_model[label].dataset["icelltype"]) > 0 + if isinstance(package, mf6.SpecificStorage) or isinstance( + package, mf6.StorageCoefficient + ): + # broadcast (layered) constants + sto_celtype = (mask * gwf_model[label].dataset["convertible"]) > 0 + if not (npf_celtype & sto_celtype).any(): + raise ValueError("Celtype need to be equal for both NPF and STO package ") + + # Check is fixed_cell option is defined for rch-package + if "fixed_cell" not in gwf_model[self.mf6_recharge_package].dataset.var(): + raise ValueError( + f"Option 'fixed_cell' is obligatory for {self.mf6_recharge_package} package, " + "when using the 'modflow_newton_formulation' of MetaMod driver" + ) + # Check if well-nodes are non-convertible + if self.mf6_wel_package is not None: + ilayer = gwf_model[self.mf6_wel_package].dataset["layer"] - 1 + irow = gwf_model[self.mf6_wel_package].dataset["row"] - 1 + icolumn = gwf_model[self.mf6_wel_package].dataset["column"] - 1 + for label in gwf_model: + package = gwf_model[label] + if isinstance(package, mf6.NodePropertyFlow): + iceltype = mask * package.dataset["icelltype"] + convertible = (iceltype[ilayer, irow, icolumn] != 0).any() + if convertible: + raise ValueError( + "Found convertible cells with irrigation extraction assigned to them" + ) def _check_sprinkling( self, msw_model: MetaSwapModel, gwf_model: GroundwaterFlowModel @@ -110,3 +156,10 @@ def write_exchanges(self, directory: Path, coupled_model: Any) -> dict[str, Any] ) return coupling_dict + + +def get_idomain(gwf_model: GroundwaterFlowModel) -> xr.DataArray: + for label in gwf_model: + package = gwf_model[label] + if isinstance(package, mf6.StructuredDiscretization): + return gwf_model[label].dataset["idomain"] diff --git a/tests/fixtures/fixture_modflow.py b/tests/fixtures/fixture_modflow.py index d2e65b78..5c5be30a 100644 --- a/tests/fixtures/fixture_modflow.py +++ b/tests/fixtures/fixture_modflow.py @@ -144,9 +144,7 @@ def make_coupled_mf6_model_newton(idomain: xr.DataArray) -> mf6.Modflow6Simulati return simulation -def make_coupled_mf6_model( - idomain: xr.DataArray, newton: bool = False -) -> mf6.Modflow6Simulation: +def make_coupled_mf6_model(idomain: xr.DataArray) -> mf6.Modflow6Simulation: _, nrow, ncol = idomain.shape gwf_model = make_mf6_model(idomain) times = get_times() From 6d4720861151e0f18abf04d67d147324185e8078 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 3 May 2024 10:48:52 +0200 Subject: [PATCH 4/7] update solver settings --- tests/fixtures/fixture_modflow.py | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/tests/fixtures/fixture_modflow.py b/tests/fixtures/fixture_modflow.py index 5c5be30a..2c9352d2 100644 --- a/tests/fixtures/fixture_modflow.py +++ b/tests/fixtures/fixture_modflow.py @@ -127,20 +127,12 @@ def make_coupled_mf6_model_newton(idomain: xr.DataArray) -> mf6.Modflow6Simulati simulation = make_mf6_simulation(gwf_model) simulation.pop("solver") - simulation["solver"] = mf6.Solution( - ["GWF_1"], - complexity="complex", - outer_dvclose=1e-5, - outer_maximum=500, - backtracking_number=0, - inner_maximum=100, - inner_dvclose=1e-6, - inner_rclose=1e-6, - rclose_option="strict", - linear_acceleration="bicgstab", - number_orthogonalizations=10, - relaxation_factor=0.0, - ) + simulation["solver"] = mf6.SolutionPresetComplex(["GWF_1"]) + simulation["solver"].dataset["outer_dvclose"] = 1e-6 + simulation["solver"].dataset["inner_dvclose"] = 1e-7 + simulation["solver"].dataset["outer_maximum"] = 500 + # don't use backtracking since we update parameters in outer iteration + simulation["solver"].dataset["backtracking_number"] = 0 return simulation From b8cee4efc6abcc9f97821f40aa4764f6f88bd586 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Wed, 12 Jun 2024 16:32:06 +0200 Subject: [PATCH 5/7] remove exchange of ss + add extra conversion to user nodes in phreatic_layer_user_indices of PhreaticBCArray --- imod_coupler/drivers/metamod/metamod.py | 29 ++++++++----------- .../kernelwrappers/mf6_newton_wrapper.py | 5 ++-- 2 files changed, 15 insertions(+), 19 deletions(-) diff --git a/imod_coupler/drivers/metamod/metamod.py b/imod_coupler/drivers/metamod/metamod.py index 56bb2a3a..233ca332 100644 --- a/imod_coupler/drivers/metamod/metamod.py +++ b/imod_coupler/drivers/metamod/metamod.py @@ -360,17 +360,18 @@ def couple(self) -> None: ) self.map_msw2mod["storage"] = conversion_matrix * self.map_msw2mod["storage"] # Create extra mapping for SY, since the SS mapping could contains a different conversion term - self.map_msw2mod["storage_sy"], self.mask_msw2mod["storage_sy"] = ( - create_mapping( - self.msw_idx, - self.node_idx, - self.msw_storage.size, - first_layer_nodes.size, - "sum", - ) + ( + self.map_msw2mod["storage_sy"], + self.mask_msw2mod["storage_sy"], + ) = create_mapping( + self.msw_idx, + self.node_idx, + self.msw_storage.size, + first_layer_nodes.size, + "sum", ) - # For exchange to SY, only the multiplication by area needs to be undone. - conversion_matrix = self.sto_conversion_terms(False, first_layer_nodes) + # For exchange to SY, act as is sc1 + conversion_matrix = self.sto_conversion_terms(True, first_layer_nodes) self.map_msw2mod["storage_sy"] = ( conversion_matrix * self.map_msw2mod["storage_sy"] ) @@ -382,7 +383,6 @@ def couple_phreatic(self, first_layer_nodes: NDArray[Any]) -> None: saturation = self.mf6.get_saturation(self.coupling.mf6_model) self.sy = self.mf6.get_sy(self.coupling.mf6_model) self.sy_top = self.sy[first_layer_nodes] - self.ss_top = self.mf6_ss[first_layer_nodes] self.recharge_nodelist = self.mf6.get_recharge_nodes( self.coupling.mf6_model, self.coupling.mf6_msw_recharge_pkg ) @@ -415,12 +415,8 @@ def exchange_msw2mod(self) -> None: self.mask_msw2mod["storage_sy"][:] * self.sy_top[:] + self.map_msw2mod["storage_sy"].dot(self.msw_storage)[:] ) - new_ss = ( - self.mask_msw2mod["storage"][:] * self.ss_top[:] - + self.map_msw2mod["storage"].dot(self.msw_storage)[:] - ) self.sto.reset() - self.sto.set(new_sy, new_ss) + self.sto.set(new_sy) self.exchange_logger.log_exchange("mf6_sy", new_sy, self.get_current_time()) self.exchange_logger.log_exchange( "msw_storage", self.msw_storage, self.get_current_time() @@ -475,7 +471,6 @@ def sto_conversion_terms( conversion_terms = 1.0 / ( self.mf6_area[indices] * (self.mf6_top[indices] - self.mf6_bot[indices]) ) - conversion_matrix = dia_matrix( (conversion_terms, [0]), shape=(self.mf6_area[indices].size, self.mf6_area[indices].size), diff --git a/imod_coupler/kernelwrappers/mf6_newton_wrapper.py b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py index d626a23a..c87d63ac 100644 --- a/imod_coupler/kernelwrappers/mf6_newton_wrapper.py +++ b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py @@ -209,6 +209,7 @@ def __init__( ptr_package_nodelist, ) self.nodes_ptr = ptr_package_nodelist # used for nodes after prepare_timestep + self.initialised = False def _set_user_indices(self) -> None: # the nodelist for bc-packages is only filled after the first prepare-timestep call @@ -235,7 +236,7 @@ def phreatic_layer_user_indices(self) -> Any: self._set_user_indices() # use initial_nodes since self.nodes is updated by set_ptr method return np.argmax(self.saturation.broadcast() > 0, axis=0).flatten()[ - self.initial_nodes + self.variable.userid[self.initial_nodes] ] @@ -284,7 +285,7 @@ def reset(self) -> None: self.sy.variable.variable_model[:] = self.initial_sy[:] self.ss.variable.variable_model[:] = self.initial_ss[:] - def set(self, new_sy: NDArray[Any], new_ss: NDArray[Any]) -> None: + def set(self, new_sy: NDArray[Any]) -> None: self.sy.set_ptr(new_sy) self.ss.set_ptr(self.zeros) From 1880209e47453b3664e5a4f50d210417825293ab Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Mon, 17 Jun 2024 11:43:05 +0200 Subject: [PATCH 6/7] change update of nodelist to after linear solve --- imod_coupler/drivers/metamod/metamod.py | 14 +++++++++++++- imod_coupler/kernelwrappers/mf6_newton_wrapper.py | 5 +++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/imod_coupler/drivers/metamod/metamod.py b/imod_coupler/drivers/metamod/metamod.py index 233ca332..a8071d18 100644 --- a/imod_coupler/drivers/metamod/metamod.py +++ b/imod_coupler/drivers/metamod/metamod.py @@ -422,7 +422,7 @@ def exchange_msw2mod(self) -> None: "msw_storage", self.msw_storage, self.get_current_time() ) - # Set recharge + # Set recharge values new_recharge = ( self.mask_msw2mod["recharge"][:] * self.mf6_recharge[:] + self.map_msw2mod["recharge"].dot(self.msw_volume)[:] / self.delt @@ -443,6 +443,18 @@ def exchange_mod2msw(self) -> None: + self.map_mod2msw["head"].dot(mf6_head)[:] ) + def do_iter(self, sol_id: int) -> bool: + """Execute a single iteration""" + self.msw.prepare_solve(0) + self.msw.solve(0) + self.exchange_msw2mod() + has_converged = self.mf6.solve(sol_id) + self.exchange_mod2msw() + self.msw.finalize_solve(0) + # update nodelist rch-package + self.rch.set_nodes() + return has_converged + def get_first_layer_nodes(self) -> NDArray[Any]: _, nrow, ncol = self.mf6.get_dis_shape(self.coupling.mf6_model) userid = self.mf6_get_userid() diff --git a/imod_coupler/kernelwrappers/mf6_newton_wrapper.py b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py index c87d63ac..de7027a3 100644 --- a/imod_coupler/kernelwrappers/mf6_newton_wrapper.py +++ b/imod_coupler/kernelwrappers/mf6_newton_wrapper.py @@ -229,6 +229,8 @@ def _set_user_indices(self) -> None: def set_ptr(self, new_values: NDArray[Any]) -> None: # The package nodelist is set to the phreatic nodes self.variable.variable_model[:] = new_values[:] + + def set_ptr_nodes(self) -> None: self.variable.nodelist[:] = (self.phreatic_modelid + 1)[:] # type ignore @property @@ -314,6 +316,9 @@ def __init__( def set(self, new_recharge: NDArray[Any]) -> None: self.recharge.set_ptr(new_recharge) + def set_nodes(self) -> None: + self.recharge.set_ptr_nodes() + class PhreaticHeads: """ From 4064b38bd9f46c7a4f06d69b8a57f90430f4497a Mon Sep 17 00:00:00 2001 From: Hendrik Kok Date: Fri, 20 Sep 2024 13:07:05 +0200 Subject: [PATCH 7/7] remove exception for non convergence --- imod_coupler/drivers/metamod/metamod.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/imod_coupler/drivers/metamod/metamod.py b/imod_coupler/drivers/metamod/metamod.py index a8071d18..ce20cfdc 100644 --- a/imod_coupler/drivers/metamod/metamod.py +++ b/imod_coupler/drivers/metamod/metamod.py @@ -253,14 +253,15 @@ def update(self) -> None: if has_converged: logger.debug(f"MF6-MSW converged in {kiter} iterations") break - if not has_converged: - logger.debug("MF6-MSW did not converge") - raise ValueError("MF6-MSW did not converge") self.mf6.finalize_solve(1) self.mf6.finalize_time_step() self.msw.finalize_time_step() + if not has_converged: + logger.debug("MF6-MSW did not converge") + # raise ValueError("MF6-MSW did not converge") + def finalize(self) -> None: self.mf6.finalize() self.msw.finalize()