A RetroSearch Logo

Home - News ( United States | United Kingdom | Italy | Germany ) - Football scores

Search Query:

Showing content from https://deltares.github.io/imod-python/user-guide/08-regridding.html below:

Regridding — iMOD Python 1.0.0rc4 documentation

Regridding#

Most MF6 packages have spatial data arrays as input. These arrays are discrete: they are defined over the simulation grid and contain values associated to each cell.

Regridding these package means: create a new package with spatial data arrays defined over a different grid. Computing what the values in these new arrays should be, is done by xugrid.

xugrid will compute a regridded array based on:

More information on the available regridding methods can be found in the xugrid documentation https://deltares.github.io/xugrid/user_guide.html

The regridding method that should be used depends on the property being regridded. For example a thermodynamically intensive property (whose value do not depend intrinsically on the grid block size) such as temperature or density can be regridded by an averaging approach (for upscaling) or sampling (for downscaling). Extensive properties (whose values do depend on the grid block size) include (water) mass and pore volume of a gridblock, and the regridding method should be chosen to reflect that. Finally regridding methods for conductivity-like properties follow the rules for parallel or serial resistors- at least when the tensor rotation angles are constant or comparable in the involved gridblocks.

Note that the different regridding methods may have a different output domain when regridding: if the original array has no-data values in some cells, then the output array may have no-data values as well, and where these end up depends on the chosen regridding method. Also note that regridding is only possible in the xy-plane, and not across the layer dimension. The output array will have the same number of layers as the input array.

Obtaining the final (i)domain#

In many real-world models, some cells will be inactive or marked as “vertical passthrough” (VPT) in the idomain array of the simulation. Some packages require that all cells that are inactictive or VPT in idomain are excluded from the package as well. An example is the npf package: cells that are inactive or VPT in idomain, should not have conductivity data in the npf package. Therefore at the end of the regridding process, a final step consists in enforcing consistency between those of idomain and all the packages. This is a 2-step process:

  1. for cells that do not have inputs in crucial packages like npf or storage, idomain will be set to inactive.

  2. for cells that are marked as inactive or VPT in idomain, all package inputs will be removed from all the packages

This synchronization step between idomain and the packages is automated, and it is carried out when calling regrid_like on the simulation object or the model object. There are 2 caveats:

  1. the synchronization between idomain and the package domains is done on the model-level. If we have a simulation containing both a flow model and a transport model then idomain for flow is determined independent from that for transport. These models may therefore end up using different domains (this may lead to undesired results, so a manual synchronization may be necessary between the flow regridding and the transport regridding) This manual synchronization can be done using the “mask_all_packages” function- this function removes input from all packages that are marked as inactive or VPT in the idomain passed to this method.

  2. The idomain/packages synchronization step is carried out when regridding a model, and when regridding a model it will use default methods for all the packages. So if you have regridded some packages yourself with non-default methods, then these are not taken into account during this synchonization step.

Regridding using default methods#

The regrid_like function is available on packages, models and simulations. When the default methods are acceptable, regridding the whole simulation is the most convenient from a user-perspective.

Regridding using non-default methods#

When non-default methods are used for one or more packages, these should be regridded separately. In that case, the most convenient approach is likely:

In code, consider an example where we want to regrid the recharge package using non default methods then we would do the following. First we’ll load some example simulation. There is a separate example contained in Regional model that you should look at if you are interested in the model building

import imod

tmpdir = imod.util.temporary_directory()

original_simulation = imod.data.hondsrug_simulation(tmpdir / "hondsrug_saved")

To reduce computational overhead for this example, we are going to clip off most of the timesteps.

original_simulation = original_simulation.clip_box(time_max="2010-01-01")

Let’s take a look at the original discretization:

original_simulation["GWF"]["dis"]

StructuredDiscretization

<xarray.Dataset> Size: 11MB
Dimensions:  (x: 500, y: 200, layer: 13)
Coordinates:
  * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 8B ...
    dy       float64 8B ...
  * layer    (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
Data variables:
    idomain  (layer, y, x) int32 5MB ...
    top      (y, x) float32 400kB ...
    bottom   (layer, y, x) float32 5MB ...

We want to regrid this to the following target grid:

import xarray as xr

import imod

dx = 100
dy = -100
xmin = 237500.0
xmax = 250000.0
ymin = 559000.0
ymax = 564000.0
target_grid = imod.util.empty_2d(
    dx=dx, dy=dy, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax
)

target_grid

<xarray.DataArray (y: 50, x: 125)> Size: 50kB
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], shape=(50, 125))
Coordinates:
  * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
  * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
    dx       float64 8B 100.0
    dy       float64 8B -100.0

This is a grid of nans, and we require a grid of ones, which can create with xarray:

target_grid = xr.ones_like(target_grid)

target_grid

<xarray.DataArray (y: 50, x: 125)> Size: 50kB
array([[1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       ...,
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.]], shape=(50, 125))
Coordinates:
  * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
  * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
    dx       float64 8B 100.0
    dy       float64 8B -100.0

Next, we’ll remove the recharge package and obtain it as a variable:

original_rch_package = original_simulation["GWF"].pop("rch")

original_rch_package

Recharge

<xarray.Dataset> Size: 21MB
Dimensions:        (y: 200, x: 500, layer: 13, time: 2)
Coordinates:
  * y              (y) float64 2kB 5.64e+05 5.64e+05 ... 5.59e+05 5.59e+05
  * x              (x) float64 4kB 2.375e+05 2.375e+05 ... 2.5e+05 2.5e+05
    dx             float64 8B ...
    dy             float64 8B ...
  * layer          (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
  * time           (time) datetime64[ns] 16B 2009-12-30T23:59:59 2009-12-31
Data variables:
    rate           (time, layer, y, x) float64 21MB nan nan nan ... nan nan nan
    print_input    bool 1B False
    print_flows    bool 1B False
    save_flows     bool 1B False
    observations   object 8B None
    repeat_stress  object 8B None

Now regrid the simulation (without recharge):

regridded_simulation = original_simulation.regrid_like("regridded", target_grid)

Let’s look at the discretization again:

regridded_simulation["GWF"]["dis"]

StructuredDiscretization

<xarray.Dataset> Size: 1MB
Dimensions:  (layer: 13, y: 50, x: 125)
Coordinates:
    dx       float64 8B 100.0
    dy       float64 8B -100.0
  * layer    (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
  * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
  * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
Data variables:
    idomain  (layer, y, x) int64 650kB -16 -16 -16 -16 -16 -16 ... 1 1 1 1 1 1
    top      (y, x) float32 25kB nan nan nan nan nan nan ... nan nan nan nan nan
    bottom   (layer, y, x) float32 325kB 6.329 6.432 6.536 ... -166.8 -166.8

All packages have been regridded, for example the NPF package:

regridded_simulation["GWF"]["npf"]

NodePropertyFlow

<xarray.Dataset> Size: 652kB
Dimensions:                              (layer: 13, y: 50, x: 125)
Coordinates:
    dx                                   float64 8B 100.0
    dy                                   float64 8B -100.0
  * layer                                (layer) int32 52B 1 2 3 4 ... 11 12 13
  * y                                    (y) float64 400B 5.64e+05 ... 5.59e+05
  * x                                    (x) float64 1kB 2.376e+05 ... 2.5e+05
Data variables: (12/22)
    icelltype                            int32 4B 0
    k                                    (layer, y, x) float32 325kB nan ... ...
    rewet                                bool 1B False
    rewet_layer                          object 8B None
    rewet_factor                         object 8B None
    rewet_iterations                     object 8B None
    ...                                   ...
    dewatered                            bool 1B True
    perched                              bool 1B True
    save_specific_discharge              bool 1B False
    save_saturation                      bool 1B False
    xt3d_option                          bool 1B False
    rhs_option                           bool 1B False

Let’s make a comparison plot of the hydraulic conductivities:

import matplotlib.pyplot as plt
import numpy as np

fig, axes = plt.subplots(nrows=2, sharex=True)
plot_kwargs = {"colors": "viridis", "levels": np.linspace(0.0, 100.0, 21), "fig": fig}

imod.visualize.spatial.plot_map(
    original_simulation["GWF"]["npf"]["k"].sel(layer=3), ax=axes[0], **plot_kwargs
)
imod.visualize.spatial.plot_map(
    regridded_simulation["GWF"]["npf"]["k"].sel(layer=3), ax=axes[1], **plot_kwargs
)

axes[0].set_ylabel("original")
axes[1].set_ylabel("regridded")

Text(28.472222222222243, 0.5, 'regridded')

Set up the input needed for custom regridding. Create a regridder weight-cache. This object can (and should) be reused for all the packages that undergo custom regridding at this stage.

from imod.util.regrid import RegridderWeightsCache

regrid_cache = RegridderWeightsCache()

regrid_cache

<imod.util.regrid.RegridderWeightsCache object at 0x0000016F386DB350>

Regrid the recharge package with a custom regridder. In this case we opt for the centroid locator regridder. This regridder is similar to using a “nearest neighbour” lookup.

from imod.common.utilities.regrid import RegridderType
from imod.mf6.regrid import RechargeRegridMethod

regridder_types = RechargeRegridMethod(rate=(RegridderType.CENTROIDLOCATOR,))

regridded_recharge = original_rch_package.regrid_like(
    target_grid,
    regrid_cache=regrid_cache,
    regridder_types=regridder_types,
)

regridded_recharge

Recharge

<xarray.Dataset> Size: 1MB
Dimensions:        (layer: 13, time: 2, y: 50, x: 125)
Coordinates:
    dx             float64 8B 100.0
    dy             float64 8B -100.0
  * layer          (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
  * time           (time) datetime64[ns] 16B 2009-12-30T23:59:59 2009-12-31
  * y              (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
  * x              (x) float64 1kB 2.376e+05 2.376e+05 ... 2.498e+05 2.5e+05
Data variables:
    rate           (time, layer, y, x) float64 1MB nan nan nan ... nan nan nan
    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

Next, add the recharge package to the regridded simulation

regridded_simulation["GWF"]["rch"] = regridded_recharge

# We can also reattach the original again
original_simulation["GWF"]["rch"] = original_rch_package

Comparison with histograms#

In the next segment we will compare the input of the models on different grids. We advice to always check how your input is regridded. In this example we upscaled grid, many input parameters are regridded with a mean method. This means that their input range is reduced, which can be seen in tailings in the histograms becoming shorter

def plot_histograms_side_by_side(array_original, array_regridded, title):
    """This function creates a plot of normalized histograms of the 2 input
    DataArray. It plots a title above each histogram."""
    _, (ax0, ax1) = plt.subplots(1, 2, sharex=True, sharey=True, tight_layout=True)
    array_original.plot.hist(ax=ax0, bins=25, density=True)
    array_regridded.plot.hist(ax=ax1, bins=25, density=True)
    ax0.title.set_text(f"{title} (original)")
    ax1.title.set_text(f"{title} (regridded)")

Compare constant head arrays.

plot_histograms_side_by_side(
    original_simulation["GWF"]["chd"]["head"],
    regridded_simulation["GWF"]["chd"]["head"],
    "chd head",
)

Compare horizontal hydraulic conductivities.

plot_histograms_side_by_side(
    original_simulation["GWF"]["npf"]["k"],
    regridded_simulation["GWF"]["npf"]["k"],
    "npf k",
)

Compare vertical hydraulic conductivities.

plot_histograms_side_by_side(
    original_simulation["GWF"]["npf"]["k33"],
    regridded_simulation["GWF"]["npf"]["k33"],
    "npf k33",
)

Compare starting heads.

plot_histograms_side_by_side(
    original_simulation["GWF"]["ic"]["start"],
    regridded_simulation["GWF"]["ic"]["start"],
    "ic start",
)

Compare river stages.

plot_histograms_side_by_side(
    original_simulation["GWF"]["riv"]["stage"],
    regridded_simulation["GWF"]["riv"]["stage"],
    "riv stage",
)

Compare river bottom elevations.

plot_histograms_side_by_side(
    original_simulation["GWF"]["riv"]["bottom_elevation"],
    regridded_simulation["GWF"]["riv"]["bottom_elevation"],
    "riv bottom elevation",
)

Compare riverbed conductance.

plot_histograms_side_by_side(
    original_simulation["GWF"]["riv"]["conductance"],
    regridded_simulation["GWF"]["riv"]["conductance"],
    "riv conductance",
)

Compare recharge rates.

plot_histograms_side_by_side(
    original_simulation["GWF"]["rch"]["rate"],
    regridded_simulation["GWF"]["rch"]["rate"],
    "rch rate",
)

Compare drainage elevations.

plot_histograms_side_by_side(
    original_simulation["GWF"]["drn-pipe"]["elevation"],
    regridded_simulation["GWF"]["drn-pipe"]["elevation"],
    "drn-pipe elevation",
)

Compare drain conductances.

plot_histograms_side_by_side(
    original_simulation["GWF"]["drn-pipe"]["conductance"],
    regridded_simulation["GWF"]["drn-pipe"]["conductance"],
    "drn-pipe conductance",
)

Let’s compare outputs now.

# Write simulation
original_simulation.write(tmpdir / "original")
# Run and open heads
original_simulation.run()
hds_original = original_simulation.open_head()

hds_original

<xarray.DataArray 'head' (time: 2, layer: 13, y: 200, x: 500)> Size: 21MB
dask.array<stack, shape=(2, 13, 200, 500), dtype=float64, chunksize=(1, 13, 200, 500), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 8B 25.0
    dy       float64 8B -25.0
  * layer    (layer) int64 104B 1 2 3 4 5 6 7 8 9 10 11 12 13
  * time     (time) float64 16B 1.157e-05 365.0

regridded_simulation.write(tmpdir / "regridded", validate=False)
regridded_simulation.run()
hds_regridded = regridded_simulation.open_head()

hds_regridded

<xarray.DataArray 'head' (time: 2, layer: 13, y: 50, x: 125)> Size: 1MB
dask.array<stack, shape=(2, 13, 50, 125), dtype=float64, chunksize=(1, 13, 50, 125), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
  * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
    dx       float64 8B 100.0
    dy       float64 8B -100.0
  * layer    (layer) int64 104B 1 2 3 4 5 6 7 8 9 10 11 12 13
  * time     (time) float64 16B 1.157e-05 365.0

Let’s make a comparison plot of the regridded heads.

fig, axes = plt.subplots(nrows=2, sharex=True)
plot_kwargs = {"colors": "viridis", "levels": np.linspace(0.0, 11.0, 12), "fig": fig}

imod.visualize.spatial.plot_map(
    hds_original.isel(layer=6, time=-1), ax=axes[0], **plot_kwargs
)
imod.visualize.spatial.plot_map(
    hds_regridded.isel(layer=6, time=-1), ax=axes[1], **plot_kwargs
)

axes[0].set_ylabel("original")
axes[1].set_ylabel("regridded")

Text(28.472222222222243, 0.5, 'regridded')

A note on regridding conductivity#

In the npf package, it is possible to use for definining the conductivity tensor:

If 1 array is given the tensor is called isotropic. Defining only K gives the same behavior as specifying K, K22 and K33 with the same value. When regridding, K33 has a default method different from that of K and K22, but it can only be applied if K33 exists in the source model in the first place. So it is recommended to introduce K33 as a separate array in the source model even if it is isotropic. Also note that default regridding methods were chosen assuming that K and K22 are roughly horizontal and K33 roughly vertical. But this may not be the case if the input arrays angle2 and/or angle3 have large values.

Regridding boundary conditions#

Special care must be taken when regridding boundary conditions, and it is recommended that users verify the balance output of a regridded simulation and compare it to the original model. If the regridded simulation is a good representation of the original simulation, the mass contributions on the balance by the different boundary conditions should be comparable in both simulations. To achieve this, it may be necessary to tweak the input or the regridding methods. An example of this is upscaling recharge (so the target grid has coarser cells than the source grid). Its default method is averaging, with the following rules:

A note on regridding transport#

Transport simulations can be unstable if constraints related to the grid Peclet number and the courant number are exceeded. This can easily happen when regridding. It may be necessary to reduce the simulation’s time step size especially when downscaling, to prevent numerical issues. Increasing dispersivities or the molecular diffusion constant can also help to stabilize the simulation. Inversely, when upscaling, a larger time step size can be acceptable.

Unsupported packages#

Some packages cannot be regridded. This includes the Lake package and the UZF package. Such packages should be removed from the simulation before regridding, and then new packages should be created by the user and then added to the regridded simulation.

This code snippet prints all default methods:

from dataclasses import asdict

import pandas as pd

from imod.tests.fixtures.package_instance_creation import ALL_PACKAGE_INSTANCES

regrid_method_setup = {
    "package name": [],
    "array name": [],
    "method name": [],
    "function name": [],
}
regrid_method_table = pd.DataFrame(regrid_method_setup)

counter = 0
for pkg in ALL_PACKAGE_INSTANCES:
    if hasattr(pkg, "_regrid_method"):
        package_name = type(pkg).__name__
        regrid_methods = asdict(pkg.get_regrid_methods())
        for array_name in regrid_methods.keys():
            method_name = regrid_methods[array_name][0].name
            function_name = ""
            if len(regrid_methods[array_name]) > 0:
                function_name = regrid_methods[array_name][1]
            regrid_method_table.loc[counter] = (
                package_name,
                array_name,
                method_name,
                function_name,
            )
            counter = counter + 1

# Set multi index to group with packages
regrid_method_table = regrid_method_table.set_index(["package name", "array name"])
# Pandas by default displays at max 60 rows, which this table exceeds.
nrows = regrid_method_table.shape[0]
# Configure pandas to increase the max amount of displayable rows.
pd.set_option("display.max_rows", nrows + 1)
# Display rows:
regrid_method_table
method name function name package name array name StructuredDiscretization top OVERLAP mean bottom OVERLAP mean idomain OVERLAP mode Dispersion diffusion_coefficient OVERLAP mean longitudinal_horizontal OVERLAP mean transversal_horizontal1 OVERLAP mean longitudinal_vertical OVERLAP mean transversal_horizontal2 OVERLAP mean transversal_vertical OVERLAP mean InitialConditions start OVERLAP mean MobileStorageTransfer porosity OVERLAP mean decay OVERLAP mean decay_sorbed OVERLAP mean bulk_density OVERLAP mean distcoef OVERLAP mean sp2 OVERLAP mean NodePropertyFlow icelltype OVERLAP mode k OVERLAP geometric_mean k22 OVERLAP geometric_mean k33 OVERLAP harmonic_mean angle1 OVERLAP mean angle2 OVERLAP mean angle3 OVERLAP mean rewet_layer OVERLAP mean SpecificStorage convertible OVERLAP mode specific_storage OVERLAP mean specific_yield OVERLAP mean StorageCoefficient convertible OVERLAP mode storage_coefficient OVERLAP mean specific_yield OVERLAP mean ConstantHead head OVERLAP mean concentration OVERLAP mean ibound OVERLAP mode Drainage elevation OVERLAP mean conductance RELATIVEOVERLAP conductance concentration OVERLAP mean Evapotranspiration surface OVERLAP mean rate OVERLAP mean depth OVERLAP mean proportion_rate OVERLAP mean proportion_depth OVERLAP mean GeneralHeadBoundary head OVERLAP mean conductance RELATIVEOVERLAP conductance concentration OVERLAP mean Recharge rate OVERLAP mean concentration OVERLAP mean River stage OVERLAP mean conductance RELATIVEOVERLAP conductance bottom_elevation OVERLAP mean concentration OVERLAP mean infiltration_factor OVERLAP mean VerticesDiscretization top OVERLAP mean bottom OVERLAP mean idomain OVERLAP mode Dispersion diffusion_coefficient OVERLAP mean longitudinal_horizontal OVERLAP mean transversal_horizontal1 OVERLAP mean longitudinal_vertical OVERLAP mean transversal_horizontal2 OVERLAP mean transversal_vertical OVERLAP mean InitialConditions start OVERLAP mean MobileStorageTransfer porosity OVERLAP mean decay OVERLAP mean decay_sorbed OVERLAP mean bulk_density OVERLAP mean distcoef OVERLAP mean sp2 OVERLAP mean NodePropertyFlow icelltype OVERLAP mode k OVERLAP geometric_mean k22 OVERLAP geometric_mean k33 OVERLAP harmonic_mean angle1 OVERLAP mean angle2 OVERLAP mean angle3 OVERLAP mean rewet_layer OVERLAP mean SpecificStorage convertible OVERLAP mode specific_storage OVERLAP mean specific_yield OVERLAP mean StorageCoefficient convertible OVERLAP mode storage_coefficient OVERLAP mean specific_yield OVERLAP mean ConstantHead head OVERLAP mean concentration OVERLAP mean ibound OVERLAP mode Drainage elevation OVERLAP mean conductance RELATIVEOVERLAP conductance concentration OVERLAP mean Evapotranspiration surface OVERLAP mean rate OVERLAP mean depth OVERLAP mean proportion_rate OVERLAP mean proportion_depth OVERLAP mean GeneralHeadBoundary head OVERLAP mean conductance RELATIVEOVERLAP conductance concentration OVERLAP mean Recharge rate OVERLAP mean concentration OVERLAP mean River stage OVERLAP mean conductance RELATIVEOVERLAP conductance bottom_elevation OVERLAP mean concentration OVERLAP mean infiltration_factor OVERLAP mean

Total running time of the script: (1 minutes 36.332 seconds)

Gallery generated by Sphinx-Gallery

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