Skip to content

Commit bae4ca3

Browse files
Issue #1632 add ATS pkg (#1645)
Fixes #1632 # Description Add ATS package. I found out that MODFLOW6 doesn't automatically forward fill stress period data for the UTL-ATS package. I guess because it is an utility? I therefore implemented forward filling, as I deem that expected behavior. Changes the following things: - Add ATS package - Add ``ats_percel`` setting to advection packages - Add ``_get_pkgkey`` method to ``Modflow6Simulation``, as packages can also be directly assigned to simulations (namely: TDIS, ATS). - Add methods to TDIS to assign and clear ATS package paths to and from its dataset. These are called upon ``Modflow6Simulation.write`` - Update outdated keys for ``inner_csvfile`` and ``outer_csvfile`` options in the imod.mf6.Solution template. - Refactor: Remove "scheme" argument from __init__ public api for advection packages, as it was only used internally (when calling regrid_like), instead make scheme a private class variable. # Checklist <!--- Before requesting review, please go through this checklist: --> - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [x] **If feature added**: Added/extended example - [x] **If feature added**: Added feature to API documentation - [ ] **If pixi.lock was changed**: Ran `pixi run generate-sbom` and committed changes
1 parent 610a3e7 commit bae4ca3

22 files changed

Lines changed: 866 additions & 71 deletions

docs/api/changelog.rst

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,12 @@ Added
2424
- Added :meth:`imod.mf6.Modflow6Simulation.create_partition_labels` to create
2525
partition labels for a MODFLOW 6 simulation from its idomain. This is useful
2626
for splitting a simulation into multiple submodels.
27+
- :class:`imod.mf6.AdaptiveTimeStepping` to specify adaptive time stepping
28+
settings for MODFLOW 6 simulations.
29+
- The ``ats_percel`` argument to :class:`imod.mf6.AdvectionTVD`,
30+
:class:`imod.mf6.AdvectionUpstream`, :class:`imod.mf6.AdvectionCentral` to
31+
adapt the time step based on the maximum fraction of a cell that a solute
32+
parcel is allowed to travel.
2733

2834
Fixed
2935
~~~~~
@@ -49,6 +55,8 @@ Fixed
4955
- Fixed bug where :meth:`imod.mf6.SourceSinkMixing.from_flow_model` would return
5056
an error upon adding a package which cannot have a ``concentration``, such as
5157
:class:`imod.mf6.HorizontalFlowBarrierResistance`.
58+
- Broken names for ``outer_csvfile`` and ``inner_csvfile`` in the
59+
:class:`imod.mf6.Solution` MODFLOW 6 template file.
5260

5361
Changed
5462
~~~~~~~

docs/api/mf6.rst

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,16 @@ Discretization
106106
TimeDiscretization.regrid_like
107107
TimeDiscretization.is_empty
108108
TimeDiscretization.get_regrid_methods
109+
AdaptiveTimeStepping
110+
AdaptiveTimeStepping.write
111+
AdaptiveTimeStepping.from_file
112+
AdaptiveTimeStepping.to_netcdf
113+
AdaptiveTimeStepping.copy
114+
AdaptiveTimeStepping.clip_box
115+
AdaptiveTimeStepping.mask
116+
AdaptiveTimeStepping.regrid_like
117+
AdaptiveTimeStepping.is_empty
118+
AdaptiveTimeStepping.get_regrid_methods
109119

110120
Model settings
111121
--------------

examples/mf6/circle_transport.py

Lines changed: 30 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
88
In overview, we'll set the following steps:
99
10-
* Setting up the flow model, just like in the circle.py example
10+
* Setting up the flow model, just like in the :doc:`circle` example
1111
* set up the transport model
1212
* Run the simulation.
1313
* Visualize the results.
@@ -134,7 +134,7 @@
134134
# Add flow model to simulation
135135
# ----------------------------
136136
#
137-
# See the :doc:`circle.py` example for more information.
137+
# See the :doc:`circle` example for more information.
138138

139139
gwf_model = imod.mf6.GroundwaterFlowModel()
140140
gwf_model["disv"] = imod.mf6.VerticesDiscretization(
@@ -184,14 +184,37 @@
184184
)
185185

186186
# %%
187+
#
188+
# Time discretization
189+
# -------------------
190+
#
187191
# Set the timesteps, we want output each year, so we specify stress periods
188-
# which last 1 year. However, timesteps of 1 year yield unstable results, so we
189-
# set ``n_timesteps`` to 10, which sets the amount of timesteps within a stress
190-
# period.
192+
# which last 1 year. We add the additional times to ensure that there is output
193+
# at the end of each year. To achieve this we setup stress periods were each
194+
# stress period is a year long. We can then use the "last" keyword in the output
195+
# control to save the output.
191196

192197
simtimes = pd.date_range(start="2000-01-01", end="2030-01-01", freq="As")
193198
simulation.create_time_discretization(additional_times=simtimes)
194-
simulation["time_discretization"]["n_timesteps"] = 10
199+
200+
# %%
201+
#
202+
# We want to use adaptive time stepping to ensure stable results without taking
203+
# an unnecessary amount of small timesteps. We set the initial time step to 0.1
204+
# day, the minimum time step to 0.1 day, and the maximum time step to 50 days.
205+
# The multiplier is set to 2.0 so that consecutive timesteps will be increased
206+
# by a factor of 2 if the convergence is good. Furthermore, we will set the
207+
# ``ats_percel`` to 0.95 in the Advection package in the next section to let
208+
# MODFLOW 6 compute an appropriate time step based on the velocity field: A
209+
# solute parcel should not travel more than one cell in a time step. From these
210+
# two criteria MODFLOW 6 will select the smallest time step.
211+
212+
simulation["ats"] = imod.mf6.AdaptiveTimeStepping(
213+
dt_init=1e-1,
214+
dt_min=1e-1,
215+
dt_max=50.0,
216+
dt_multiplier=2.0,
217+
)
195218

196219
# %%
197220
# Buoyancy
@@ -241,7 +264,7 @@
241264
xt3d_off=False,
242265
xt3d_rhs=False,
243266
)
244-
transport_model["adv"] = imod.mf6.AdvectionUpstream()
267+
transport_model["adv"] = imod.mf6.AdvectionTVD(ats_percel=0.95)
245268
transport_model["mst"] = imod.mf6.MobileStorageTransfer(porosity)
246269

247270
# %%

examples/user-guide/07-time-discretization.py

Lines changed: 63 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -32,8 +32,10 @@
3232
import imod
3333

3434
# %%
35+
#
3536
# We can discretize a simulation as follows, this creates a
36-
# TimeDiscretization object under the key ``"time_discretization"``.
37+
# :class:`imod.mf6.TimeDiscretization` object under the key
38+
# ``"time_discretization"``.
3739

3840
simulation = imod.mf6.Modflow6Simulation("example")
3941
simulation.create_time_discretization(
@@ -43,7 +45,8 @@
4345
simulation["time_discretization"]
4446

4547
# %%
46-
# To view the data inside TimeDiscretization object:
48+
#
49+
# To view the data inside :class:`imod.mf6.TimeDiscretization` object:
4750

4851
simulation["time_discretization"].dataset
4952

@@ -128,15 +131,15 @@
128131
# Notice that iMOD Python figured out that the two boundary conditions, both
129132
# with two timesteps, should lead to three stress periods!
130133
#
131-
# Specifying extra settings
132-
# -------------------------
134+
# Extra settings
135+
# --------------
133136
#
134-
# The ``TimeDiscretization`` package also supports other settings, like the
135-
# amount of timesteps which Modflow should use within a stress period, as well
136-
# as a timestep multiplier, to gradually increase the timesteps modflow uses
137-
# within a stress period. This can be useful when boundary conditions change
138-
# very abruptly between stress periods. These settings are set by accessing
139-
# their respective variables in the ``dataset``.
137+
# The :class:`imod.mf6.TimeDiscretization` package also supports other settings,
138+
# like the amount of timesteps which Modflow should use within a stress period,
139+
# as well as a timestep multiplier, to gradually increase the timesteps modflow
140+
# uses within a stress period. This can be useful when boundary conditions
141+
# change very abruptly between stress periods. These settings are set by
142+
# accessing their respective variables in the ``dataset``.
140143

141144
times = simulation_bc["time_discretization"].dataset.coords["time"]
142145

@@ -148,3 +151,53 @@
148151
simulation_bc["time_discretization"].dataset
149152

150153
# %%
154+
#
155+
# Adaptive Time Stepping
156+
# ----------------------
157+
#
158+
# The :class:`imod.mf6.AdaptiveTimeStepping` package can be used to let Modflow
159+
# decide the length of each timestep within a stress period. This can be useful
160+
# when boundary conditions change abruptly, or when a model is highly dynamic.
161+
# The :class:`imod.mf6.AdaptiveTimeStepping` package can be added to a
162+
# :class:`imod.mf6.Modflow6Simulation` as follows:
163+
164+
simulation_bc["ats"] = imod.mf6.AdaptiveTimeStepping(
165+
dt_init=1e-4, dt_min=1e-4, dt_max=2.0, dt_multiplier=2.0, dt_fail_divisor=5.0
166+
)
167+
168+
simulation_bc["ats"]
169+
170+
# %%
171+
#
172+
# The ``dt_init`` parameter sets the initial timestep length, ``dt_min`` and
173+
# ``dt_max`` set the minimum and maximum timestep length allowed,
174+
# ``dt_multiplier`` multiplies the timestep length for each consecutive timestep
175+
# if the model converges well, and ``dt_fail_divisor`` decreases the timestep
176+
# length if the model has convergence problems.
177+
#
178+
# There is one further way to influence the timestep length, which is to set the
179+
# ``ats_percel`` parameter in the :class:`imod.mf6.AdvectionTVD` package of a
180+
# transport model. This parameter sets the maximum fraction of a cell that a
181+
# parcel of solute is allowed to travel in a single timestep. This can be used
182+
# to ensure that the timestep is small enough to get accurate transport results.
183+
# Once this is set, MODFLOW 6 will choose the smallest necessary timestep based
184+
# on the convergence criteria and the ``ats_percel`` criteria.
185+
186+
transport_model = imod.mf6.GroundwaterTransportModel()
187+
transport_model["adv"] = imod.mf6.AdvectionTVD(ats_percel=0.95)
188+
simulation_bc["gwt_1"] = transport_model
189+
190+
transport_model["adv"]
191+
192+
# %%
193+
#
194+
# Conclusion
195+
# ----------
196+
#
197+
# The :meth:`imod.mf6.Modflow6Simulation.create_time_discretization` method can
198+
# be used to easily generate stress periods for a simulation, based on the
199+
# timesteps of the boundary conditions assigned to the model. The
200+
# :class:`imod.mf6.AdaptiveTimeStepping` package can be used to let MODFLOW 6
201+
# decide the length of each timestep within a stress period, based on
202+
# convergence criteria and, in case of using the ``ats_percel`` options, the
203+
# velocity field.

imod/mf6/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
from imod.mf6.adv import AdvectionCentral, AdvectionTVD, AdvectionUpstream
66
from imod.mf6.api_package import ApiPackage
7+
from imod.mf6.ats import AdaptiveTimeStepping
78
from imod.mf6.buy import Buoyancy
89
from imod.mf6.chd import ConstantHead
910
from imod.mf6.cnc import ConstantConcentration

imod/mf6/adv.py

Lines changed: 105 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -9,49 +9,97 @@
99
robust alternative.
1010
"""
1111

12+
from abc import ABC
1213
from copy import deepcopy
14+
from typing import Optional
15+
16+
import numpy as np
1317

1418
from imod.common.interfaces.iregridpackage import IRegridPackage
19+
from imod.common.utilities.dataclass_type import DataclassType
1520
from imod.mf6.package import Package
21+
from imod.schemata import AllValueSchema, DimsSchema, DTypeSchema
22+
from imod.typing import GridDataArray
23+
from imod.util.regrid import RegridderWeightsCache
1624

1725

18-
class Advection(Package, IRegridPackage):
26+
class AdvectionBase(Package, IRegridPackage, ABC):
1927
_pkg_id = "adv"
2028
_template = Package._initialize_template(_pkg_id)
21-
22-
def __init__(self, scheme: str):
23-
dict_dataset = {"scheme": scheme}
29+
_scheme: str
30+
31+
_init_schemata = {
32+
"ats_percel": [
33+
DimsSchema(),
34+
DTypeSchema(np.floating),
35+
AllValueSchema(">", 0.0),
36+
],
37+
}
38+
39+
def __init__(self, ats_percel: Optional[float] = None, validate: bool = True):
40+
dict_dataset = {"scheme": self._scheme, "ats_percel": ats_percel}
2441
super().__init__(dict_dataset)
42+
self._validate_init_schemata(validate)
2543

2644
def _render(self, directory, pkgname, globaltimes, binary):
27-
scheme = self.dataset["scheme"].item()
28-
return self._template.render({"scheme": scheme})
45+
render_dict = {}
46+
render_dict["scheme"] = self.dataset["scheme"].item()
47+
if "ats_percel" in self.dataset and self._valid(
48+
self.dataset["ats_percel"].item()
49+
):
50+
render_dict["ats_percel"] = self.dataset["ats_percel"].item()
51+
return self._template.render(render_dict)
2952

3053
def mask(self, _) -> Package:
3154
"""
32-
The mask method is irrelevant for this package , instead this method
55+
The mask method is irrelevant for this package, instead this method
56+
retuns a copy of itself.
57+
"""
58+
return deepcopy(self)
59+
60+
def regrid_like(
61+
self,
62+
target_grid: GridDataArray,
63+
regrid_cache: RegridderWeightsCache,
64+
regridder_types: Optional[DataclassType] = None,
65+
) -> Package:
66+
"""
67+
The regrid_like method is irrelevant for this package, instead this method
3368
retuns a copy of itself.
3469
"""
3570
return deepcopy(self)
3671

3772

38-
class AdvectionUpstream(Advection):
73+
class AdvectionUpstream(AdvectionBase):
3974
"""
4075
The upstream weighting (first order upwind) scheme sets the concentration
4176
at the cellface between two adjacent cells equal to the concentration in
4277
the cell where the flow comes from. It surpresses oscillations.
43-
Note: all constructor arguments will be ignored
78+
79+
Parameters
80+
----------
81+
ats_percel: float, optional
82+
Fractional cell distance submitted by the ADV Package to the
83+
:class:`imod.mf6.AdaptiveTimeStepping` (ATS) package. If ``ats_percel``
84+
is specified and the ATS Package is active, a time step calculation will
85+
be made for each cell based on flow through the cell and cell
86+
properties. The largest time step will be calculated such that the
87+
advective fractional cell distance (``ats_percel``) is not exceeded for
88+
any active cell in the grid. This time-step constraint will be submitted
89+
to the ATS Package, perhaps with constraints submitted by other
90+
packages, in the calculation of the time step. ``ats_percel`` must be
91+
greater than zero. If a value of zero is specified for ``ats_percel``
92+
the program will automatically reset it to an internal no data value to
93+
indicate that time steps should not be subject to this constraint.
94+
Requires MODFLOW 6.6.0 or higher.
95+
validate: bool, optional
96+
Validate the package upon initialization. Defaults to True.
4497
"""
4598

46-
def __init__(self, scheme: str = "upstream"):
47-
if not scheme == "upstream":
48-
raise ValueError(
49-
"error in scheme parameter. Should be 'upstream' if present."
50-
)
51-
super().__init__(scheme="upstream")
99+
_scheme = "upstream"
52100

53101

54-
class AdvectionCentral(Advection):
102+
class AdvectionCentral(AdvectionBase):
55103
"""
56104
The central-in-space weighting scheme is based on a simple
57105
distance-weighted linear interpolation between the center of cell n and the
@@ -60,25 +108,53 @@ class AdvectionCentral(Advection):
60108
grids without equal spacing between connected cells, it is retained here
61109
for consistency with nomenclature used by other MODFLOW-based transport
62110
programs, such as MT3D.
63-
Note: all constructor arguments will be ignored
111+
112+
Parameters
113+
----------
114+
ats_percel: float, optional
115+
Fractional cell distance submitted by the ADV Package to the
116+
:class:`imod.mf6.AdaptiveTimeStepping` (ATS) package. If ``ats_percel``
117+
is specified and the ATS Package is active, a time step calculation will
118+
be made for each cell based on flow through the cell and cell
119+
properties. The largest time step will be calculated such that the
120+
advective fractional cell distance (``ats_percel``) is not exceeded for
121+
any active cell in the grid. This time-step constraint will be submitted
122+
to the ATS Package, perhaps with constraints submitted by other
123+
packages, in the calculation of the time step. ``ats_percel`` must be
124+
greater than zero. If a value of zero is specified for ``ats_percel``
125+
the program will automatically reset it to an internal no data value to
126+
indicate that time steps should not be subject to this constraint.
127+
Requires MODFLOW 6.6.0 or higher.
128+
validate: bool, optional
129+
Validate the package upon initialization. Defaults to True.
64130
"""
65131

66-
def __init__(self, scheme: str = "central"):
67-
if not scheme == "central":
68-
raise ValueError(
69-
"error in scheme parameter. Should be 'central' if present."
70-
)
71-
super().__init__(scheme="central")
132+
_scheme = "central"
72133

73134

74-
class AdvectionTVD(Advection):
135+
class AdvectionTVD(AdvectionBase):
75136
"""
76137
An implicit second order TVD scheme. More expensive than upstream
77138
weighting but more robust.
78-
Note: all constructor arguments will be ignored
139+
140+
Parameters
141+
----------
142+
ats_percel: float, optional
143+
Fractional cell distance submitted by the ADV Package to the
144+
:class:`imod.mf6.AdaptiveTimeStepping` (ATS) package. If ``ats_percel``
145+
is specified and the ATS Package is active, a time step calculation will
146+
be made for each cell based on flow through the cell and cell
147+
properties. The largest time step will be calculated such that the
148+
advective fractional cell distance (``ats_percel``) is not exceeded for
149+
any active cell in the grid. This time-step constraint will be submitted
150+
to the ATS Package, perhaps with constraints submitted by other
151+
packages, in the calculation of the time step. ``ats_percel`` must be
152+
greater than zero. If a value of zero is specified for ``ats_percel``
153+
the program will automatically reset it to an internal no data value to
154+
indicate that time steps should not be subject to this constraint.
155+
Requires MODFLOW 6.6.0 or higher.
156+
validate: bool, optional
157+
Validate the package upon initialization. Defaults to True.
79158
"""
80159

81-
def __init__(self, scheme: str = "TVD"):
82-
if not scheme == "TVD":
83-
raise ValueError("error in scheme parameter. Should be 'TVD' if present.")
84-
super().__init__(scheme="TVD")
160+
_scheme = "TVD"

0 commit comments

Comments
 (0)