From aad7e89f7019d5fc17b88db7ec91c616193f1bde Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Mon, 6 May 2024 16:00:15 +0200 Subject: [PATCH 1/8] Add support for discharge scaling in drainage package --- imod/mf6/boundary_condition.py | 2 +- imod/mf6/drn.py | 42 +++++++++++++++++++++++++++++++++- imod/templates/mf6/gwf-drn.j2 | 2 ++ 3 files changed, 44 insertions(+), 2 deletions(-) diff --git a/imod/mf6/boundary_condition.py b/imod/mf6/boundary_condition.py index 8f1466a59..e5c4b6e60 100644 --- a/imod/mf6/boundary_condition.py +++ b/imod/mf6/boundary_condition.py @@ -303,7 +303,7 @@ def get_period_varnames(self): if hasattr(self, "_auxiliary_data"): result.extend(get_variable_names(self)) - return result + return [var for var in result if var in self.dataset.data_vars] class AdvancedBoundaryCondition(BoundaryCondition, abc.ABC): diff --git a/imod/mf6/drn.py b/imod/mf6/drn.py index 66a4d868e..6b6095636 100644 --- a/imod/mf6/drn.py +++ b/imod/mf6/drn.py @@ -3,6 +3,7 @@ import numpy as np from imod.logging import init_log_decorator +from imod.mf6.auxiliary_variables import get_variable_names from imod.mf6.boundary_condition import BoundaryCondition from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.mf6.utilities.regrid import RegridderType @@ -37,6 +38,15 @@ class Drainage(BoundaryCondition, IRegridPackage): concentration_boundary_type: ({"AUX", "AUXMIXED"}, optional) if this flow package is used in simulations also involving transport, then this keyword specifies how outflow over this boundary is computed. + scaling_depth: array of floats (xr.DataArray) + if provided, enables the drainage package discharge scaling. + Effectively, the conductance is reduced. Its value may be both positive + or negative. + If positive: drainage starts at head equals elevation. At (elevation + + depth) the scaling factor becomes 1.0. + If negative: drainage starts at + elevation + depth (elevation - |depth|). At elevation the scaling + factor becomes 1.0. print_input: ({True, False}, optional) keyword to indicate that the list of drain information will be written to the listing file immediately after it is read. Default is False. @@ -93,6 +103,12 @@ class Drainage(BoundaryCondition, IRegridPackage): ), CONC_DIMS_SCHEMA, ], + "scaling_depth": [ + DTypeSchema(np.floating), + IndexesSchema(), + CoordsSchema(("layer",)), + BOUNDARY_DIMS_SCHEMA, + ], "print_flows": [DTypeSchema(np.bool_), DimsSchema()], "save_flows": [DTypeSchema(np.bool_), DimsSchema()], } @@ -104,9 +120,10 @@ class Drainage(BoundaryCondition, IRegridPackage): ], "conductance": [IdentityNoDataSchema("elevation"), AllValueSchema(">", 0.0)], "concentration": [IdentityNoDataSchema("elevation"), AllValueSchema(">=", 0.0)], + "scaling_depth": [IdentityNoDataSchema("elevation")], } - _period_data = ("elevation", "conductance") + _period_data = ("elevation", "conductance", "scaling_depth") _keyword_map = {} _template = BoundaryCondition._initialize_template(_pkg_id) _auxiliary_data = {"concentration": "species"} @@ -123,6 +140,7 @@ def __init__( elevation, conductance, concentration=None, + scaling_depth=None, concentration_boundary_type="aux", print_input=False, print_flows=False, @@ -135,6 +153,7 @@ def __init__( "elevation": elevation, "conductance": conductance, "concentration": concentration, + "scaling_depth": scaling_depth, "concentration_boundary_type": concentration_boundary_type, "print_input": print_input, "print_flows": print_flows, @@ -153,5 +172,26 @@ def _validate(self, schemata, **kwargs): return errors + def render(self, directory, pkgname, globaltimes, binary): + """Render fills in the template only, doesn't write binary data""" + d = {"binary": binary} + bin_ds = self._get_bin_ds() + d["periods"] = self._period_paths( + directory, pkgname, globaltimes, bin_ds, binary + ) + # construct the rest (dict for render) + d = self._get_options(d) + d["maxbound"] = self._max_active_n() + + if (hasattr(self, "_auxiliary_data")) and (names := get_variable_names(self)): + d["auxiliary"] = names + if self.dataset["scaling_depth"] is not None: + d["auxdepthname"] = "scaling_depth" + if "auxiliary" in d: + d["auxiliary"] += ["scaling_depth"] + else: + d["auxiliary"] = ["scaling_depth"] + return self._template.render(d) + def get_regrid_methods(self) -> Optional[dict[str, Tuple[RegridderType, str]]]: return self._regrid_method diff --git a/imod/templates/mf6/gwf-drn.j2 b/imod/templates/mf6/gwf-drn.j2 index 3771dd22a..73dc5f494 100644 --- a/imod/templates/mf6/gwf-drn.j2 +++ b/imod/templates/mf6/gwf-drn.j2 @@ -3,6 +3,8 @@ begin options {% endif %} {%- if auxmultname is defined %} auxmultname {{auxmultname}} {% endif %} +{%- if auxdepthname is defined %} auxdepthname {{auxdepthname}} +{% endif %} {%- if boundnames is defined %} boundnames {% endif %} {%- if print_input is defined %} print_input From fc6fbecde538df70edd74fe5279c947dc6284654 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Tue, 7 May 2024 14:50:12 +0200 Subject: [PATCH 2/8] Update ims options MODFLOW 6 has changed a number of options in the IMS package (since quite some time). * This adds support for outer_csvfile and inner_csvfile instead of csv_output. * This simplifies the no_ptc option: instead of a bool and a possible additional entry, just the entry is required. (In the IMS file, NO_PTC is equal to NO_PCT ALL, i.e. the ALL is an assumed default value.) * This also adds support for the adaptive time stepping entry ats_outer_maximum_fraction. All the deprecated solution keyword arguments have been removed from the examples. --- examples/mf6/Henry.py | 4 - examples/mf6/circle.py | 2 - examples/mf6/circle_transport.py | 4 - examples/mf6/ex01_twri.py | 2 - examples/mf6/example_1d_transport.py | 4 - examples/mf6/example_models.py | 4 - examples/mf6/hondsrug.py | 2 - examples/mf6/lake.py | 2 - examples/mf6/transport_2d.py | 4 - imod/mf6/ims.py | 87 +++++++++++-------- imod/templates/mf6/sln-ims.j2 | 9 +- .../flow_transport_simulation_fixture.py | 4 - imod/tests/fixtures/mf6_circle_fixture.py | 8 -- imod/tests/fixtures/mf6_modelrun_fixture.py | 2 - .../fixtures/mf6_rectangle_with_lakes.py | 2 - .../fixtures/mf6_small_models_fixture.py | 2 - imod/tests/fixtures/mf6_twri_disv_fixture.py | 2 - imod/tests/fixtures/mf6_twri_fixture.py | 2 - imod/tests/fixtures/msw_model_fixture.py | 2 - imod/tests/test_benchmark.py | 3 +- .../test_formats/test_disv_conversion.py | 2 - imod/tests/test_mf6/test_mf6_ims.py | 51 +++++++++-- imod/tests/test_mf6/test_mf6_uzf_model.py | 2 - 23 files changed, 106 insertions(+), 100 deletions(-) diff --git a/examples/mf6/Henry.py b/examples/mf6/Henry.py index 745d090e2..2e3c0d8e9 100644 --- a/examples/mf6/Henry.py +++ b/examples/mf6/Henry.py @@ -205,8 +205,6 @@ simulation["flow_solver"] = imod.mf6.Solution( modelnames=["flow"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-6, outer_maximum=500, under_relaxation=None, @@ -221,8 +219,6 @@ simulation["transport_solver"] = imod.mf6.Solution( modelnames=["transport"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-6, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/circle.py b/examples/mf6/circle.py index 8bdedb3d1..3bb534097 100644 --- a/examples/mf6/circle.py +++ b/examples/mf6/circle.py @@ -136,8 +136,6 @@ simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/circle_transport.py b/examples/mf6/circle_transport.py index 427d9f97a..75acf2045 100644 --- a/examples/mf6/circle_transport.py +++ b/examples/mf6/circle_transport.py @@ -175,8 +175,6 @@ simulation["flow_solver"] = imod.mf6.Solution( modelnames=["flow"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -264,8 +262,6 @@ simulation["transport_solver"] = imod.mf6.Solution( modelnames=["transport"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/ex01_twri.py b/examples/mf6/ex01_twri.py index 5f47908c8..e690701a0 100644 --- a/examples/mf6/ex01_twri.py +++ b/examples/mf6/ex01_twri.py @@ -161,8 +161,6 @@ simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/example_1d_transport.py b/examples/mf6/example_1d_transport.py index 760d6c327..0576e7a28 100644 --- a/examples/mf6/example_1d_transport.py +++ b/examples/mf6/example_1d_transport.py @@ -171,8 +171,6 @@ def create_transport_model(flowmodel, speciesname, dispersivity, retardation, de simulation["flow_solver"] = imod.mf6.Solution( modelnames=["flow"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -187,8 +185,6 @@ def create_transport_model(flowmodel, speciesname, dispersivity, retardation, de simulation["transport_solver"] = imod.mf6.Solution( modelnames=["tpt_a", "tpt_b", "tpt_c", "tpt_d"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/example_models.py b/examples/mf6/example_models.py index 35001ada8..b9f39e420 100644 --- a/examples/mf6/example_models.py +++ b/examples/mf6/example_models.py @@ -193,8 +193,6 @@ def create_twri_simulation() -> imod.mf6.Modflow6Simulation: simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -281,8 +279,6 @@ def create_circle_simulation(): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/hondsrug.py b/examples/mf6/hondsrug.py index dd71508f0..d9a21c20b 100644 --- a/examples/mf6/hondsrug.py +++ b/examples/mf6/hondsrug.py @@ -577,8 +577,6 @@ def outer_edge(da): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/lake.py b/examples/mf6/lake.py index efdd40c10..c2572e438 100644 --- a/examples/mf6/lake.py +++ b/examples/mf6/lake.py @@ -164,8 +164,6 @@ simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/examples/mf6/transport_2d.py b/examples/mf6/transport_2d.py index 6d14be7e7..2700773ef 100644 --- a/examples/mf6/transport_2d.py +++ b/examples/mf6/transport_2d.py @@ -150,8 +150,6 @@ simulation["flow_solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -166,8 +164,6 @@ simulation["transport_solver"] = imod.mf6.Solution( modelnames=["TPT_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/mf6/ims.py b/imod/mf6/ims.py index 904b59070..1c4d7e017 100644 --- a/imod/mf6/ims.py +++ b/imod/mf6/ims.py @@ -285,9 +285,6 @@ class Solution(Package): scaling method in Hill (1992). l2norm - symmetric matrix scaling using the L2 norm. Default value: None - SolutionPresetSimple: None - SolutionPresetModerate: None - SolutionPresetComplex: None reordering_method: str options: {"None", "rcm", "md"} an optional keyword that defines the matrix reordering approach used. By @@ -296,9 +293,6 @@ class Solution(Package): rcm - reverse Cuthill McKee ordering. md - minimum degree ordering Default value: None - SolutionPresetSimple: None - SolutionPresetModerate: None - SolutionPresetComplex: None print_option: str options: {"None", "summary", "all"} is a flag that controls printing of convergence information from the @@ -311,17 +305,15 @@ class Solution(Package): convergence information to each model listing file in addition to SUMMARY information. Default value: "summary" - SolutionPresetSimple: No Default - SolutionPresetModerate: No Default - SolutionPresetComplex: No Default - csv_output: str, optional - False if no csv is to be written for the output, enter str of filename + outer_csvfile: str, optional + None if no csv is to be written for the output, str of filename + if csv is to be written. + Default value: None + inner_csvfile: str, optional + None if no csv is to be written for the output, str of filename if csv is to be written. - Default value: False - SolutionPresetSimple: No Default - SolutionPresetModerate: No Default - SolutionPresetComplex: No Default - no_ptc: ({True, False}, optional) + Default value: None + no_ptc: str, optional, either None, "all", or "first". is a flag that is used to disable pseudo-transient continuation (PTC). Option only applies to steady-state stress periods for models using the Newton-Raphson formulation. For many problems, PTC can significantly @@ -333,10 +325,22 @@ class Solution(Package): PTC option should be used to deactivate PTC. This NO PTC option should also be used in order to compare convergence behavior with other MODFLOW versions, as PTC is only available in MODFLOW 6. - Default value: False - SolutionPresetSimple: No Default - SolutionPresetModerate: No Default - SolutionPresetComplex: No Default + ats_outer_maximum_fraction: float, optional. + real value defining the fraction of the maximum allowable outer iterations + used with the Adaptive Time Step (ATS) capability if it is active. If this + value is set to zero by the user, then this solution will have no effect + on ATS behavior. This value must be greater than or equal to zero and less + than or equal to 0.5 or the program will terminate with an error. If it is + not specified by the user, then it is assigned a default value of one + third. When the number of outer iterations for this solution is less than + the product of this value and the maximum allowable outer iterations, then + ATS will increase the time step length by a factor of DTADJ in the ATS + input file. When the number of outer iterations for this solution is + greater than the maximum allowable outer iterations minus the product of + this value and the maximum allowable outer iterations, then the ATS (if + active) will decrease the time step length by a factor of 1 / DTADJ. + + Default value is None. validate: {True, False} Flag to indicate whether the package should be validated upon initialization. This raises a ValidationError if package input is @@ -391,8 +395,10 @@ def __init__( scaling_method=None, reordering_method=None, print_option="summary", - csv_output=False, - no_ptc=False, + outer_csvfile=None, + inner_csvfile=None, + no_ptc=None, + ats_outer_maximum_fraction=None, validate: bool = True, ): dict_dataset = { @@ -419,7 +425,9 @@ def __init__( "scaling_method": scaling_method, "reordering_method": reordering_method, "print_option": print_option, - "csv_output": csv_output, + "outer_csvfile": outer_csvfile, + "inner_csvfile": inner_csvfile, + "ats_outer_maximum_fraction": ats_outer_maximum_fraction, "no_ptc": no_ptc, } # Make sure the modelnames are set as a variable rather than dimension: @@ -460,12 +468,17 @@ def get_models_in_solution(self) -> list[str]: def SolutionPresetSimple( - modelnames, print_option="summary", csv_output=False, no_ptc=False + modelnames, + print_option="summary", + outer_csvfile=None, + inner_csvfile=None, + no_ptc=None, ): solution = Solution( modelnames=modelnames, print_option=print_option, - csv_output=csv_output, + outer_csvfile=outer_csvfile, + inner_csvfile=inner_csvfile, no_ptc=no_ptc, outer_dvclose=0.001, outer_maximum=25, @@ -487,19 +500,22 @@ def SolutionPresetSimple( preconditioner_levels=0, preconditioner_drop_tolerance=0, number_orthogonalizations=0, - scaling_method=None, - reordering_method=None, ) return solution def SolutionPresetModerate( - modelnames, print_option="summary", csv_output=False, no_ptc=False + modelnames, + print_option="summary", + outer_csvfile=None, + inner_csvfile=None, + no_ptc=None, ): solution = Solution( modelnames=modelnames, print_option=print_option, - csv_output=csv_output, + outer_csvfile=outer_csvfile, + inner_csvfile=inner_csvfile, no_ptc=no_ptc, outer_dvclose=0.01, outer_maximum=50, @@ -521,19 +537,22 @@ def SolutionPresetModerate( preconditioner_levels=0, preconditioner_drop_tolerance=0.0, number_orthogonalizations=0, - scaling_method=None, - reordering_method=None, ) return solution def SolutionPresetComplex( - modelnames, print_option="summary", csv_output=False, no_ptc=False + modelnames, + print_option="summary", + outer_csvfile=None, + inner_csvfile=None, + no_ptc=None, ): solution = Solution( modelnames=modelnames, print_option=print_option, - csv_output=csv_output, + outer_csvfile=outer_csvfile, + inner_csvfile=inner_csvfile, no_ptc=no_ptc, outer_dvclose=0.1, outer_maximum=100, @@ -555,7 +574,5 @@ def SolutionPresetComplex( preconditioner_levels=5, preconditioner_drop_tolerance=0.0001, number_orthogonalizations=2, - scaling_method=None, - reordering_method=None, ) return solution diff --git a/imod/templates/mf6/sln-ims.j2 b/imod/templates/mf6/sln-ims.j2 index f8ad6ce0d..52531ff26 100644 --- a/imod/templates/mf6/sln-ims.j2 +++ b/imod/templates/mf6/sln-ims.j2 @@ -1,10 +1,15 @@ begin options {% if print_option is defined %} print_option {{print_option}} {% endif %} -{%- if csv_output_filerecord is defined %} csv_output fileout {{csvfile}} +{%- if outer_csvfile is defined %} outer_csvfile fileout {{outer_csvfile}} {% endif %} -{%- if no_ptcrecord is defined %} no_ptc {% if no_ptc_option is defined %}{{no_ptc_option}}{% endif %} +{%- if inner_csvfile is defined %} inner_csvfile fileout {{inner_csvfile}} +{% endif %} +{%- if no_ptc is defined %} no_ptc {{no_ptc}} +{% endif -%} +{%- if ats_outer_maximum_fraction is defined %} ats_outer_maximum_fraction {{ats_outer_maximum_fraction}} {% endif -%} + end options begin nonlinear diff --git a/imod/tests/fixtures/flow_transport_simulation_fixture.py b/imod/tests/fixtures/flow_transport_simulation_fixture.py index b8f6fff4e..67f54d37b 100644 --- a/imod/tests/fixtures/flow_transport_simulation_fixture.py +++ b/imod/tests/fixtures/flow_transport_simulation_fixture.py @@ -181,8 +181,6 @@ def flow_transport_simulation(): simulation["solver"] = imod.mf6.Solution( modelnames=["flow"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -197,8 +195,6 @@ def flow_transport_simulation(): simulation["transport_solver"] = imod.mf6.Solution( modelnames=["tpt_a", "tpt_b", "tpt_c", "tpt_d"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-6, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_circle_fixture.py b/imod/tests/fixtures/mf6_circle_fixture.py index 6e7b6a89b..2c4c38b37 100644 --- a/imod/tests/fixtures/mf6_circle_fixture.py +++ b/imod/tests/fixtures/mf6_circle_fixture.py @@ -63,8 +63,6 @@ def make_circle_model(): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -148,8 +146,6 @@ def make_circle_model_flow_with_transport_data(species: list[str]): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -329,8 +325,6 @@ def circle_model_transport(): simulation["transport_solver"] = imod.mf6.Solution( modelnames=["transport"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -408,8 +402,6 @@ def circle_model_transport_multispecies_variable_density(): simulation["transport_solver"] = imod.mf6.Solution( modelnames=modelnames, print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_modelrun_fixture.py b/imod/tests/fixtures/mf6_modelrun_fixture.py index acaced6b4..1605bbc1c 100644 --- a/imod/tests/fixtures/mf6_modelrun_fixture.py +++ b/imod/tests/fixtures/mf6_modelrun_fixture.py @@ -58,8 +58,6 @@ def assert_model_can_run( simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_rectangle_with_lakes.py b/imod/tests/fixtures/mf6_rectangle_with_lakes.py index abc8d593a..fbc27be99 100644 --- a/imod/tests/fixtures/mf6_rectangle_with_lakes.py +++ b/imod/tests/fixtures/mf6_rectangle_with_lakes.py @@ -78,8 +78,6 @@ def rectangle_with_lakes(): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_small_models_fixture.py b/imod/tests/fixtures/mf6_small_models_fixture.py index fb8587a94..00bffc4b2 100644 --- a/imod/tests/fixtures/mf6_small_models_fixture.py +++ b/imod/tests/fixtures/mf6_small_models_fixture.py @@ -164,8 +164,6 @@ def solution_settings() -> imod.mf6.Solution: return imod.mf6.Solution( modelnames=["flow"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_twri_disv_fixture.py b/imod/tests/fixtures/mf6_twri_disv_fixture.py index 7dad3c2b6..4170aaa20 100644 --- a/imod/tests/fixtures/mf6_twri_disv_fixture.py +++ b/imod/tests/fixtures/mf6_twri_disv_fixture.py @@ -136,8 +136,6 @@ def toface(da): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/mf6_twri_fixture.py b/imod/tests/fixtures/mf6_twri_fixture.py index 92adb35b0..c884ca8fe 100644 --- a/imod/tests/fixtures/mf6_twri_fixture.py +++ b/imod/tests/fixtures/mf6_twri_fixture.py @@ -191,8 +191,6 @@ def make_twri_model(): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/fixtures/msw_model_fixture.py b/imod/tests/fixtures/msw_model_fixture.py index 943a90fb6..3d16262af 100644 --- a/imod/tests/fixtures/msw_model_fixture.py +++ b/imod/tests/fixtures/msw_model_fixture.py @@ -91,8 +91,6 @@ def make_coupled_mf6_model(): simulation["solver"] = mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/test_benchmark.py b/imod/tests/test_benchmark.py index cfc378266..a77e7d3ef 100644 --- a/imod/tests/test_benchmark.py +++ b/imod/tests/test_benchmark.py @@ -165,7 +165,8 @@ def setup_mf6_basic_simulation_imod(): simulation["gwf"] = gwf_model # Define solver settings simulation["solver"] = imod.mf6.SolutionPresetSimple( - modelnames=["gwf"], print_option="summary", csv_output=False, no_ptc=True + modelnames=["gwf"], + print_option="summary", ) # Collect time discretization diff --git a/imod/tests/test_formats/test_disv_conversion.py b/imod/tests/test_formats/test_disv_conversion.py index 647c220d0..48e11fabc 100644 --- a/imod/tests/test_formats/test_disv_conversion.py +++ b/imod/tests/test_formats/test_disv_conversion.py @@ -49,8 +49,6 @@ def test_convert_to_disv(imodflow_model, tmp_path, create_grid): simulation["solver"] = imod.mf6.Solution( modelnames=["gwf"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, diff --git a/imod/tests/test_mf6/test_mf6_ims.py b/imod/tests/test_mf6/test_mf6_ims.py index 9bdb40b38..031135a6e 100644 --- a/imod/tests/test_mf6/test_mf6_ims.py +++ b/imod/tests/test_mf6/test_mf6_ims.py @@ -10,8 +10,6 @@ def create_ims() -> imod.mf6.Solution: return imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -29,8 +27,6 @@ def test_render(): ims = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, @@ -71,8 +67,6 @@ def test_wrong_dtype(): imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=4, outer_maximum=500, under_relaxation=None, @@ -105,3 +99,48 @@ def test_add_already_present_model(): ims = create_ims() with pytest.raises(ValueError): ims.add_model_to_solution("GWF_1") + + +def test_ims_options(): + """ + Ensure inner/outer csvfile, no_ptc, and ats_outer_maximum_fraction are + written. + """ + ims = imod.mf6.Solution( + modelnames=["GWF_1"], + outer_dvclose=0.001, + outer_maximum=20, + inner_maximum=100, + inner_dvclose=0.0001, + inner_rclose=1.0, + linear_acceleration="cg", + inner_csvfile="inner.csv", + outer_csvfile="outer.csv", + no_ptc="first", + ats_outer_maximum_fraction=0.25, + ) + actual = ims.render(None, None, None, None) + expected = textwrap.dedent( + """\ + begin options + print_option summary + outer_csvfile fileout outer.csv + inner_csvfile fileout inner.csv + no_ptc first + ats_outer_maximum_fraction 0.25 + end options + + begin nonlinear + outer_dvclose 0.001 + outer_maximum 20 + end nonlinear + + begin linear + inner_maximum 100 + inner_dvclose 0.0001 + inner_rclose 1.0 + linear_acceleration cg + end linear + """ + ) + assert expected == actual diff --git a/imod/tests/test_mf6/test_mf6_uzf_model.py b/imod/tests/test_mf6/test_mf6_uzf_model.py index 557761383..a34843027 100644 --- a/imod/tests/test_mf6/test_mf6_uzf_model.py +++ b/imod/tests/test_mf6/test_mf6_uzf_model.py @@ -135,8 +135,6 @@ def uzf_model(): simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", - csv_output=False, - no_ptc=True, outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, From 5174f2a05af8ac197ca09306f36944de616b6f07 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Tue, 7 May 2024 14:52:00 +0200 Subject: [PATCH 3/8] Formatting --- imod/mf6/ims.py | 2 +- imod/tests/test_mf6/test_mf6_ims.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/mf6/ims.py b/imod/mf6/ims.py index 1c4d7e017..ad527fa38 100644 --- a/imod/mf6/ims.py +++ b/imod/mf6/ims.py @@ -338,7 +338,7 @@ class Solution(Package): input file. When the number of outer iterations for this solution is greater than the maximum allowable outer iterations minus the product of this value and the maximum allowable outer iterations, then the ATS (if - active) will decrease the time step length by a factor of 1 / DTADJ. + active) will decrease the time step length by a factor of 1 / DTADJ. Default value is None. validate: {True, False} diff --git a/imod/tests/test_mf6/test_mf6_ims.py b/imod/tests/test_mf6/test_mf6_ims.py index 031135a6e..8f79089ce 100644 --- a/imod/tests/test_mf6/test_mf6_ims.py +++ b/imod/tests/test_mf6/test_mf6_ims.py @@ -106,7 +106,7 @@ def test_ims_options(): Ensure inner/outer csvfile, no_ptc, and ats_outer_maximum_fraction are written. """ - ims = imod.mf6.Solution( + ims = imod.mf6.Solution( modelnames=["GWF_1"], outer_dvclose=0.001, outer_maximum=20, From b99db7a4e773a037563d4d63cf60228e0405a3c0 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Wed, 8 May 2024 11:18:12 +0200 Subject: [PATCH 4/8] Add OptionSchema and apply it for mf6 Solution --- imod/mf6/ims.py | 32 +++++++++++++++-------- imod/schemata.py | 26 ++++++++++++++++++- imod/tests/test_mf6/test_mf6_ims.py | 38 ++++++++++++++++++++++++++++ imod/tests/test_mf6/test_schemata.py | 30 ++++++++++++++++++++++ 4 files changed, 115 insertions(+), 11 deletions(-) create mode 100644 imod/tests/test_mf6/test_schemata.py diff --git a/imod/mf6/ims.py b/imod/mf6/ims.py index ad527fa38..88a2221e8 100644 --- a/imod/mf6/ims.py +++ b/imod/mf6/ims.py @@ -3,7 +3,7 @@ from imod.logging import init_log_decorator from imod.mf6.package import Package -from imod.schemata import DTypeSchema +from imod.schemata import AllValueSchema, DTypeSchema, OptionSchema class Solution(Package): @@ -36,7 +36,7 @@ class Solution(Package): SolutionPresetComplex: 0.1 outer_maximum: int integer value defining the maximum number of outer (nonlinear) - iterations – that is, calls to the solution routine. For a linear + iterations - that is, calls to the solution routine. For a linear problem outer_maximum should be 1. SolutionPresetSimple: 25 SolutionPresetModerate: 50 @@ -81,7 +81,7 @@ class Solution(Package): SolutionPresetModerate: "bicgstab" SolutionPresetComplex: "bicgstab" under_relaxation: str, optional - options: {"None", "simple", "cooley", "bdb"} + options: {None, "simple", "cooley", "bdb"} is an optional keyword that defines the nonlinear relative_rclose schemes used. By default under_relaxation is not used. None - relative_rclose is not used. @@ -178,7 +178,7 @@ class Solution(Package): backtracking_reduction_factor: float, optional real value defining the reduction in step size used for residual reduction computations. The value of backtracking reduction factor is - between 142 MODFLOW 6 – Description of Input and Output zero and one. + between 142 MODFLOW 6 - Description of Input and Output zero and one. The value usually ranges from 0.1 to 0.3; a value of 0.2 works well for most problems. backtracking_reduction_factor only needs to be specified if backtracking number is greater than zero. @@ -202,17 +202,17 @@ class Solution(Package): options: {"strict", "l2norm_rclose", "relative_rclose"} an optional keyword that defines the specific flow residual criterion used. - strict– an optional keyword that is used to specify that inner rclose + strict: an optional keyword that is used to specify that inner rclose represents a infinity-norm (absolute convergence criteria) and that the head and flow convergence criteria must be met on the first inner iteration (this criteria is equivalent to the criteria used by the MODFLOW-2005 PCG package (Hill, 1990)). - l2norm_rclose – an optionalkeyword that is used to specify that inner + l2norm_rclose: an optionalkeyword that is used to specify that inner rclose represents a l-2 norm closure criteria instead of a infinity-norm (absolute convergence criteria). When l2norm_rclose is specified, a reasonable initial inner rclose value is 0.1 times the number of active cells when meters and seconds are the defined MODFLOW 6 length and time. - relative_rclose – an optional keyword that is used to specify that + relative_rclose: an optional keyword that is used to specify that inner_rclose represents a relative L-2 Norm reduction closure criteria instead of a infinity-Norm (absolute convergence criteria). When relative_rclose is specified, a reasonable initial inner_rclose value is @@ -277,7 +277,7 @@ class Solution(Package): SolutionPresetModerate: 0 SolutionPresetComplex: 2 scaling_method: str - options: {"None", "diagonal", "l2norm"} + options: {None, "diagonal", "l2norm"} an optional keyword that defines the matrix scaling approach used. By default, matrix scaling is not applied. None - no matrix scaling applied. @@ -286,7 +286,7 @@ class Solution(Package): l2norm - symmetric matrix scaling using the L2 norm. Default value: None reordering_method: str - options: {"None", "rcm", "md"} + options: {None, "rcm", "md"} an optional keyword that defines the matrix reordering approach used. By default, matrix reordering is not applied. None - original ordering. @@ -294,7 +294,7 @@ class Solution(Package): md - minimum degree ordering Default value: None print_option: str - options: {"None", "summary", "all"} + options: {None, "summary", "all"} is a flag that controls printing of convergence information from the solver. None - means print nothing. @@ -356,6 +356,9 @@ class Solution(Package): "inner_maximum": [DTypeSchema(np.integer)], "inner_dvclose": [DTypeSchema(np.floating)], "inner_rclose": [DTypeSchema(np.floating)], + "linear_acceleration": [OptionSchema(("cg", "bicgstab"))], + "rclose_option": [OptionSchema(("strict", "l2norm_rclose", "relative_rclose"))], + "under_relaxation": [OptionSchema(("simple", "cooley", "bdb"))], "under_relaxation_theta": [DTypeSchema(np.floating)], "under_relaxation_kappa": [DTypeSchema(np.floating)], "under_relaxation_gamma": [DTypeSchema(np.floating)], @@ -365,6 +368,15 @@ class Solution(Package): "backtracking_reduction_factor": [DTypeSchema(np.floating)], "backtracking_residual_limit": [DTypeSchema(np.floating)], "number_orthogonalizations": [DTypeSchema(np.integer)], + "scaling_method": [OptionSchema(("diagonal", "l2norm"))], + "reordering_method": [OptionSchema(("rcm", "md"))], + "print_option": [OptionSchema(("summary", "all"))], + "no_ptc": [OptionSchema(("first", "all"))], + "ats_outer_maximum_fraction": [ + DTypeSchema(np.floating), + AllValueSchema(">=", 0.0), + AllValueSchema("<=", 0.5), + ], } _template = Package._initialize_template(_pkg_id) diff --git a/imod/schemata.py b/imod/schemata.py index 5aae2b157..e71600be6 100644 --- a/imod/schemata.py +++ b/imod/schemata.py @@ -34,7 +34,7 @@ import abc import operator from functools import partial -from typing import Any, Callable, Dict, Optional, Tuple, TypeAlias, Union +from typing import Any, Callable, Dict, Optional, Sequence, Tuple, TypeAlias, Union import numpy as np import scipy @@ -148,6 +148,30 @@ def __or__(self, other): return SchemaUnion(*self.schemata, other) +class OptionSchema(BaseSchema): + """ + Check whether the value is one of given valid options. + """ + + def __init__(self, options: Sequence[Any]): + self.options = options + + def validate(self, obj: ScalarAsDataArray, **kwargs) -> None: + if scalar_None(obj): + return + + # MODFLOW 6 is not case sensitive for string options. + value = obj.item() + if isinstance(value, str): + value = value.lower() + + if value not in self.options: + valid_values = ", ".join(map(str, self.options)) + raise ValidationError( + f"Invalid option: {value}. Valid options are: {valid_values}" + ) + + class DTypeSchema(BaseSchema): def __init__(self, dtype: DTypeLike) -> None: if dtype in [ diff --git a/imod/tests/test_mf6/test_mf6_ims.py b/imod/tests/test_mf6/test_mf6_ims.py index 8f79089ce..7f4397fff 100644 --- a/imod/tests/test_mf6/test_mf6_ims.py +++ b/imod/tests/test_mf6/test_mf6_ims.py @@ -1,3 +1,4 @@ +import re import textwrap import pytest @@ -144,3 +145,40 @@ def test_ims_options(): """ ) assert expected == actual + + +def test_ims_option_validation(): + expected = textwrap.dedent( + """ + * linear_acceleration + \t- Invalid option: abc. Valid options are: cg, bicgstab + * rclose_option + \t- Invalid option: any. Valid options are: strict, l2norm_rclose, relative_rclose + * scaling_method + \t- Invalid option: random. Valid options are: diagonal, l2norm + * reordering_method + \t- Invalid option: alphabetical. Valid options are: rcm, md + * print_option + \t- Invalid option: whatever. Valid options are: summary, all + * no_ptc + \t- Invalid option: last. Valid options are: first, all + * ats_outer_maximum_fraction + \t- not all values comply with criterion: <= 0.5""" + ) + + with pytest.raises(ValidationError, match=re.escape(expected)): + imod.mf6.Solution( + modelnames=["GWF_1"], + outer_dvclose=0.001, + outer_maximum=20, + inner_maximum=100, + inner_dvclose=0.0001, + inner_rclose=1.0, + rclose_option="any", + linear_acceleration="abc", + scaling_method="random", + reordering_method="alphabetical", + print_option="whatever", + no_ptc="last", + ats_outer_maximum_fraction=1.0, + ) diff --git a/imod/tests/test_mf6/test_schemata.py b/imod/tests/test_mf6/test_schemata.py new file mode 100644 index 000000000..7b321d584 --- /dev/null +++ b/imod/tests/test_mf6/test_schemata.py @@ -0,0 +1,30 @@ +import pytest +import xarray as xr + +from imod import schemata as sch +from imod.schemata import ValidationError + + +def test_option_schema(): + # Test for integer values + schema = sch.OptionSchema([1, 2]) + # None values should be skipped + assert schema.validate(xr.DataArray(None)) is None + assert schema.validate(xr.DataArray(1)) is None + assert schema.validate(xr.DataArray(2)) is None + + with pytest.raises(ValidationError): + schema.validate(xr.DataArray(3)) + + # Test for string values + schema = sch.OptionSchema(["a", "b"]) + # None values should be skipped + assert schema.validate(xr.DataArray(None)) is None + assert schema.validate(xr.DataArray("a")) is None + # String values are treated case insensitive + assert schema.validate(xr.DataArray("B")) is None + + with pytest.raises( + ValidationError, match="Invalid option: c. Valid options are: a, b" + ): + schema.validate(xr.DataArray("c")) From 0f696662c4d729786f99655a7d413102c4912b8c Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Wed, 8 May 2024 11:37:32 +0200 Subject: [PATCH 5/8] Allow print_option "none" string. Remove some non-ASCII chars from docstring --- imod/mf6/ims.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/mf6/ims.py b/imod/mf6/ims.py index 88a2221e8..955840826 100644 --- a/imod/mf6/ims.py +++ b/imod/mf6/ims.py @@ -216,7 +216,7 @@ class Solution(Package): inner_rclose represents a relative L-2 Norm reduction closure criteria instead of a infinity-Norm (absolute convergence criteria). When relative_rclose is specified, a reasonable initial inner_rclose value is - 1.0 × 10−4 and convergence is achieved for a given inner (linear) + 1.0 * 10-4 and convergence is achieved for a given inner (linear) iteration when ∆h ≤ inner_dvclose and the current L-2 Norm is ≤ the product of the relativ_rclose and the initial L-2 Norm for the current inner (linear) iteration. If rclose_option is not specified, an absolute @@ -294,7 +294,7 @@ class Solution(Package): md - minimum degree ordering Default value: None print_option: str - options: {None, "summary", "all"} + options: {"none", "summary", "all"} is a flag that controls printing of convergence information from the solver. None - means print nothing. @@ -370,7 +370,7 @@ class Solution(Package): "number_orthogonalizations": [DTypeSchema(np.integer)], "scaling_method": [OptionSchema(("diagonal", "l2norm"))], "reordering_method": [OptionSchema(("rcm", "md"))], - "print_option": [OptionSchema(("summary", "all"))], + "print_option": [OptionSchema(("none", "summary", "all"))], "no_ptc": [OptionSchema(("first", "all"))], "ats_outer_maximum_fraction": [ DTypeSchema(np.floating), From e917981efca3b91467e0a06a6fc438e09041e3ed Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Wed, 8 May 2024 11:42:29 +0200 Subject: [PATCH 6/8] Update changelog --- docs/api/changelog.rst | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index e79114efb..c6bc4cd7d 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -29,7 +29,10 @@ Added for phreatic models is to use the Newton formulation, where cells remain active, and this option irrelevant. Option Added in :func:`imod.mf6.recharge`, - +- Added support for ats_outer_maximum_fraction in :class:`imod.mf6.Solution`. +- Added validation for linear_acceleration, rclose_option, scaling_method, + reordering_method, print_option and no_ptc entries in + :class:`imod.mf6.Solution`. Fixed ~~~~~ @@ -41,6 +44,12 @@ Fixed ipf and idf files are now applied - Fix bug where y-coords were flipped in :class:`imod.msw.MeteoMapping` +Changed +~~~~~~~ +- Replaced csv_output by outer_csvfile and inner_csvfile in + :class:`imod.mf6.Solution` to match newer MODFLOW 6 releases. +- Changed no_ptc from a bool to an option string in :class:`imod.mf6.Solution`. + [0.16.0] - 2024-03-29 --------------------- From 56fc6db80266bffb63d12f6df24638d171671598 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Wed, 8 May 2024 12:16:32 +0200 Subject: [PATCH 7/8] bdb -> dbd ... --- imod/mf6/ims.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/mf6/ims.py b/imod/mf6/ims.py index 955840826..193809dbf 100644 --- a/imod/mf6/ims.py +++ b/imod/mf6/ims.py @@ -81,7 +81,7 @@ class Solution(Package): SolutionPresetModerate: "bicgstab" SolutionPresetComplex: "bicgstab" under_relaxation: str, optional - options: {None, "simple", "cooley", "bdb"} + options: {None, "simple", "cooley", "dbd"} is an optional keyword that defines the nonlinear relative_rclose schemes used. By default under_relaxation is not used. None - relative_rclose is not used. @@ -358,7 +358,7 @@ class Solution(Package): "inner_rclose": [DTypeSchema(np.floating)], "linear_acceleration": [OptionSchema(("cg", "bicgstab"))], "rclose_option": [OptionSchema(("strict", "l2norm_rclose", "relative_rclose"))], - "under_relaxation": [OptionSchema(("simple", "cooley", "bdb"))], + "under_relaxation": [OptionSchema(("simple", "cooley", "dbd"))], "under_relaxation_theta": [DTypeSchema(np.floating)], "under_relaxation_kappa": [DTypeSchema(np.floating)], "under_relaxation_gamma": [DTypeSchema(np.floating)], From 82a08d45d14b710bb14a3beade9ba8c5e49241e1 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Wed, 8 May 2024 13:36:28 +0200 Subject: [PATCH 8/8] Forgot to update test for print option one... --- imod/tests/test_mf6/test_mf6_ims.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/tests/test_mf6/test_mf6_ims.py b/imod/tests/test_mf6/test_mf6_ims.py index 7f4397fff..48ab9d77f 100644 --- a/imod/tests/test_mf6/test_mf6_ims.py +++ b/imod/tests/test_mf6/test_mf6_ims.py @@ -159,7 +159,7 @@ def test_ims_option_validation(): * reordering_method \t- Invalid option: alphabetical. Valid options are: rcm, md * print_option - \t- Invalid option: whatever. Valid options are: summary, all + \t- Invalid option: whatever. Valid options are: none, summary, all * no_ptc \t- Invalid option: last. Valid options are: first, all * ats_outer_maximum_fraction