iMOD Python provides nice functionality to discretize your models into stress periods, depending on the timesteps you assigned your boundary conditions. This functionality is activated with the create_time_discretization()
method.
To demonstrate the create_time_discretization()
method, we first have to create a Model object. In this case we’ll use a MODFLOW 6 simulation, but the imod.wq.SeawatModel
also supports this.
Note
The simulation created in this example misses some mandatory packes, so is not able to write a functional simulation. It is purely intended to describe the workings of the time_discretization
method. For a fully functional MODFLOW 6 model, see the examples in MODFLOW 6.
Wel’ll start off with the usual imports:
import numpy as np import pandas as pd import xarray as xr import imod
We can discretize a simulation as follows, this creates a TimeDiscretization object under the key "time_discretization"
.
simulation = imod.mf6.Modflow6Simulation("example") simulation.create_time_discretization( additional_times=["2000-01-01", "2000-01-02", "2000-01-04"] ) simulation["time_discretization"]
TimeDiscretization
<xarray.Dataset> Size: 48B Dimensions: (time: 2) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-02 Data variables: timestep_duration (time) float64 16B 1.0 2.0 n_timesteps int64 8B 1 timestep_multiplier float64 8B 1.0
time
(time)
datetime64[ns]
2000-01-01 2000-01-02
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000'], dtype='datetime64[ns]')
timestep_duration
(time)
float64
1.0 2.0
n_timesteps
()
int64
1
timestep_multiplier
()
float64
1.0
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', name='time', freq=None))
To view the data inside TimeDiscretization object:
simulation["time_discretization"].dataset
<xarray.Dataset> Size: 48B Dimensions: (time: 2) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-02 Data variables: timestep_duration (time) float64 16B 1.0 2.0 n_timesteps int64 8B 1 timestep_multiplier float64 8B 1.0
time
(time)
datetime64[ns]
2000-01-01 2000-01-02
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000'], dtype='datetime64[ns]')
timestep_duration
(time)
float64
1.0 2.0
n_timesteps
()
int64
1
timestep_multiplier
()
float64
1.0
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', name='time', freq=None))
Notice that even though we specified three points in time, only two timesteps are included in the time
coordinate, this is because Modflow requires a start time and a duration of each stress period. iMOD Python therefore uses three points in time to compute two stress periods with a duration.
simulation["time_discretization"].dataset["timestep_duration"]
<xarray.DataArray 'timestep_duration' (time: 2)> Size: 16B array([1., 2.]) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-02
time
(time)
datetime64[ns]
2000-01-01 2000-01-02
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000'], dtype='datetime64[ns]')
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', name='time', freq=None))
These two stress periods use their respective start time in their time
coordinate.
The create_time_discretization
method becomes especially useful if we add boundary conditions to our groundwater model. We’ll first still have to initialize a groundwater flow model though:
gwf_model = imod.mf6.GroundwaterFlowModel()
We’ll create a dummy grid, which is used to provide an appropriate grid to the boundary condition packages.
template_data = np.ones((2, 1, 2, 2)) layer = [1] y = [1.5, 0.5] x = [0.5, 1.5] dims = ("time", "layer", "y", "x")
Next, we can assign a Constant Head boundary with two timesteps:
chd_times = pd.to_datetime(["2000-01-01", "2000-01-04"]) chd_data = xr.DataArray( data=template_data, dims=dims, coords={"time": chd_times, "layer": layer, "y": y, "x": x}, ) gwf_model["chd"] = imod.mf6.ConstantHead(head=chd_data) gwf_model["chd"].dataset
<xarray.Dataset> Size: 139B Dimensions: (time: 2, layer: 1, y: 2, x: 2) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-04 * layer (layer) int64 8B 1 * y (y) float64 16B 1.5 0.5 * x (x) float64 16B 0.5 1.5 Data variables: head (time, layer, y, x) float64 64B 1.0 1.0 1.0 ... 1.0 1.0 1.0 print_input bool 1B False print_flows bool 1B False save_flows bool 1B False observations object 8B None repeat_stress object 8B None
time
(time)
datetime64[ns]
2000-01-01 2000-01-04
array(['2000-01-01T00:00:00.000000000', '2000-01-04T00:00:00.000000000'], dtype='datetime64[ns]')
layer
(layer)
int64
1
y
(y)
float64
1.5 0.5
x
(x)
float64
0.5 1.5
head
(time, layer, y, x)
float64
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
array([[[[1., 1.], [1., 1.]]], [[[1., 1.], [1., 1.]]]])
print_input
()
bool
False
print_flows
()
bool
False
save_flows
()
bool
False
observations
()
object
None
array(None, dtype=object)
repeat_stress
()
object
None
array(None, dtype=object)
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-04'], dtype='datetime64[ns]', name='time', freq=None))
PandasIndex
PandasIndex(Index([1], dtype='int64', name='layer'))
PandasIndex
PandasIndex(Index([1.5, 0.5], dtype='float64', name='y'))
PandasIndex
PandasIndex(Index([0.5, 1.5], dtype='float64', name='x'))
We’ll also assign a Recharge boundary with two timesteps, which differ from the ConstantHead boundary:
rch_times = pd.to_datetime(["2000-01-01", "2000-01-02"]) rch_data = xr.DataArray( data=template_data, dims=dims, coords={"time": rch_times, "layer": layer, "y": y, "x": x}, ) gwf_model["rch"] = imod.mf6.Recharge(rate=rch_data) gwf_model["rch"].dataset
<xarray.Dataset> Size: 140B Dimensions: (time: 2, layer: 1, y: 2, x: 2) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-02 * layer (layer) int64 8B 1 * y (y) float64 16B 1.5 0.5 * x (x) float64 16B 0.5 1.5 Data variables: rate (time, layer, y, x) float64 64B 1.0 1.0 1.0 ... 1.0 1.0 1.0 print_input bool 1B False print_flows bool 1B False save_flows bool 1B False observations object 8B None repeat_stress object 8B None fixed_cell bool 1B False
time
(time)
datetime64[ns]
2000-01-01 2000-01-02
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000'], dtype='datetime64[ns]')
layer
(layer)
int64
1
y
(y)
float64
1.5 0.5
x
(x)
float64
0.5 1.5
rate
(time, layer, y, x)
float64
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
array([[[[1., 1.], [1., 1.]]], [[[1., 1.], [1., 1.]]]])
print_input
()
bool
False
print_flows
()
bool
False
save_flows
()
bool
False
observations
()
object
None
array(None, dtype=object)
repeat_stress
()
object
None
array(None, dtype=object)
fixed_cell
()
bool
False
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', name='time', freq=None))
PandasIndex
PandasIndex(Index([1], dtype='int64', name='layer'))
PandasIndex
PandasIndex(Index([1.5, 0.5], dtype='float64', name='y'))
PandasIndex
PandasIndex(Index([0.5, 1.5], dtype='float64', name='x'))
We can now let iMOD Python figure out how the simulation’s time should be discretized. It is important that we provide an endtime, otherwise the duration of the last stress period cannot be determined:
endtime = pd.to_datetime(["2000-01-06"]) simulation_bc = imod.mf6.Modflow6Simulation("example_bc") simulation_bc["gwf_1"] = gwf_model simulation_bc.create_time_discretization(additional_times=endtime) simulation_bc["time_discretization"].dataset
<xarray.Dataset> Size: 64B Dimensions: (time: 3) Coordinates: * time (time) datetime64[ns] 24B 2000-01-01 ... 2000-01-04 Data variables: timestep_duration (time) float64 24B 1.0 2.0 2.0 n_timesteps int64 8B 1 timestep_multiplier float64 8B 1.0
time
(time)
datetime64[ns]
2000-01-01 2000-01-02 2000-01-04
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000', '2000-01-04T00:00:00.000000000'], dtype='datetime64[ns]')
timestep_duration
(time)
float64
1.0 2.0 2.0
n_timesteps
()
int64
1
timestep_multiplier
()
float64
1.0
PandasIndex
PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02', '2000-01-04'], dtype='datetime64[ns]', name='time', freq=None))
Notice that iMOD Python figured out that the two boundary conditions, both with two timesteps, should lead to three stress periods!
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4