.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\circle.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_mf6_circle.py: Circle ====== This example illustrates how to setup a very simple unstructured groundwater model using the ``imod`` package and associated packages. In overview, we'll set the following steps: * Create a triangular mesh for a disk geometry. * Create the xugrid UgridDataArrays containg the MODFLOW6 parameters. * Feed these arrays into the imod mf6 classes. * Write to modflow6 files. * Run the model. * Open the results back into UgridDataArrays. * Visualize the results. .. GENERATED FROM PYTHON SOURCE LINES 18-20 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 22-23 We'll start with the following imports: .. GENERATED FROM PYTHON SOURCE LINES 23-32 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import xarray as xr import xugrid as xu from pandas import isnull import imod .. GENERATED FROM PYTHON SOURCE LINES 33-41 Create a mesh ------------- The first steps consists of generating a mesh. In this example, we'll use data included with iMOD Python for a circular mesh. Note that this is a `Ugrid2D object. `_ For more information on working with unstructured grids see the `Xugrid documentation `_ .. GENERATED FROM PYTHON SOURCE LINES 41-47 .. code-block:: Python grid = imod.data.circle() grid .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 10kB Dimensions: (mesh2d_nFaces: 216, mesh2d_nMax_face_nodes: 3, mesh2d_nEdges: 342, two: 2, mesh2d_nNodes: 127) Coordinates: mesh2d_node_x (mesh2d_nNodes) float64 1kB 0.0 6.123e-14 ... 331.1 157.5 mesh2d_node_y (mesh2d_nNodes) float64 1kB 0.0 1e+03 ... 764.7 818.3 Dimensions without coordinates: mesh2d_nFaces, mesh2d_nMax_face_nodes, mesh2d_nEdges, two, mesh2d_nNodes Data variables: mesh2d int64 8B 0 mesh2d_face_nodes (mesh2d_nFaces, mesh2d_nMax_face_nodes) int32 3kB 0 ... 1 mesh2d_edge_nodes (mesh2d_nEdges, two) int64 5kB 0 7 0 17 ... 126 125 126 Attributes: Conventions: CF-1.9 UGRID-1.0 .. GENERATED FROM PYTHON SOURCE LINES 48-49 We can plot this object as follows: .. GENERATED FROM PYTHON SOURCE LINES 49-54 .. code-block:: Python fig, ax = plt.subplots() xu.plot.line(grid, ax=ax) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_001.png :alt: circle :srcset: /examples/mf6/images/sphx_glr_circle_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 55-71 Create UgridDataArray --------------------- Now that we have defined the grid, we can start defining the model parameter data. Our goal here is to define a steady-state model with: * Uniform conductivities of 1.0 m/d; * Two layers of 5.0 m thick; * Uniform recharge of 0.001 m/d on the top layer; * Constant heads of 1.0 m along the exterior edges of the mesh. From these boundary conditions, we would expect circular mounding of the groundwater; with small flows in the center and larger flows as the recharge accumulates while the groundwater flows towards the exterior boundary. .. GENERATED FROM PYTHON SOURCE LINES 71-89 .. code-block:: Python nface = grid.n_face nlayer = 2 idomain = xu.UgridDataArray( xr.DataArray( np.ones((nlayer, nface), dtype=np.int32), coords={"layer": [1, 2]}, dims=["layer", grid.face_dimension], ), grid=grid, ) icelltype = xu.full_like(idomain, 0) k = xu.full_like(idomain, 1.0, dtype=float) k33 = k.copy() rch_rate = xu.full_like(idomain.sel(layer=1), 0.001, dtype=float) bottom = idomain * xr.DataArray([5.0, 0.0], dims=["layer"]) .. GENERATED FROM PYTHON SOURCE LINES 90-93 All the data above have been constants over the grid. For the constant head boundary, we'd like to only set values on the external border. We can `py:method:xugrid.UgridDataset.binary_dilation` to easily find these cells: .. GENERATED FROM PYTHON SOURCE LINES 93-104 .. code-block:: Python chd_location = xu.zeros_like(idomain.sel(layer=2), dtype=bool).ugrid.binary_dilation( border_value=True ) constant_head = xu.full_like(idomain.sel(layer=2), 1.0, dtype=float).where(chd_location) fig, ax = plt.subplots() constant_head.ugrid.plot(ax=ax) xu.plot.line(grid, ax=ax, color="black") ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_002.png :alt: layer = 2 :srcset: /examples/mf6/images/sphx_glr_circle_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 105-110 Write the model --------------- The first step is to define an empty model, the parameters and boundary conditions are added in the form of the familiar MODFLOW packages. .. GENERATED FROM PYTHON SOURCE LINES 110-152 .. code-block:: Python gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["disv"] = imod.mf6.VerticesDiscretization( top=10.0, bottom=bottom, idomain=idomain ) gwf_model["chd"] = imod.mf6.ConstantHead( constant_head, print_input=True, print_flows=True, save_flows=True ) gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0) gwf_model["npf"] = imod.mf6.NodePropertyFlow( icelltype=icelltype, k=k, k33=k33, save_flows=True, ) gwf_model["sto"] = imod.mf6.SpecificStorage( specific_storage=1.0e-5, specific_yield=0.15, transient=False, convertible=0, ) gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all") gwf_model["rch"] = imod.mf6.Recharge(rch_rate) simulation = imod.mf6.Modflow6Simulation("circle") simulation["GWF_1"] = gwf_model simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, inner_dvclose=1.0e-4, inner_rclose=0.001, inner_maximum=100, linear_acceleration="cg", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"]) .. GENERATED FROM PYTHON SOURCE LINES 153-154 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 154-158 .. code-block:: Python modeldir = imod.util.temporary_directory() simulation.write(modeldir) .. GENERATED FROM PYTHON SOURCE LINES 159-167 Run the model ------------- .. note:: The following lines assume the ``mf6`` executable is available on your PATH. :ref:`The Modflow 6 examples introduction ` shortly describes how to add it to yours. .. GENERATED FROM PYTHON SOURCE LINES 167-170 .. code-block:: Python simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 171-175 Open the results ---------------- First, we'll open the heads (.hds) file. .. GENERATED FROM PYTHON SOURCE LINES 175-180 .. code-block:: Python head = simulation.open_head() head .. raw:: html 
<xarray.DataArray 'head' (time: 1, layer: 2, mesh2d_nFaces: 216)> Size: 3kB
    dask.array<stack, shape=(1, 2, 216), dtype=float64, chunksize=(1, 2, 216), chunktype=numpy.ndarray>
    Coordinates:
      * layer          (layer) int64 16B 1 2
      * time           (time) float64 8B 1.0
      * mesh2d_nFaces  (mesh2d_nFaces) int64 2kB 0 1 2 3 4 5 ... 211 212 213 214 215


.. GENERATED FROM PYTHON SOURCE LINES 181-186 For a DISV MODFLOW6 model, the heads are returned as a UgridDataArray. While all layers are timesteps are available, they are only loaded into memory as needed. We may also open the cell-by-cell flows (.cbc) file. .. GENERATED FROM PYTHON SOURCE LINES 186-191 .. code-block:: Python cbc = simulation.open_flow_budget() print(cbc.keys()) .. rst-class:: sphx-glr-script-out .. code-block:: none KeysView( Size: 28kB Dimensions: (layer: 2, time: 1, mesh2d_nEdges: 342, mesh2d_nFaces: 216) Coordinates: * layer (layer) int64 16B 1 2 * time (time) float64 8B 1.0 * mesh2d_nEdges (mesh2d_nEdges) int64 3kB 0 1 2 3 ... 339 340 341 * mesh2d_nFaces (mesh2d_nFaces) int64 2kB 0 1 2 3 ... 213 214 215 Data variables: flow-horizontal-face (time, layer, mesh2d_nEdges) float64 5kB dask.array flow-horizontal-face-x (time, layer, mesh2d_nEdges) float64 5kB dask.array flow-horizontal-face-y (time, layer, mesh2d_nEdges) float64 5kB dask.array flow-lower-face (time, layer, mesh2d_nFaces) float64 3kB dask.array chd_chd (time, layer, mesh2d_nFaces) float64 3kB dask.array) .. GENERATED FROM PYTHON SOURCE LINES 192-201 The flows are returned as a dictionary of UgridDataArrays. This dictionary contains all entries that are stored in the CBC file, but like for the heads file the data are only loaded into memory when needed. The horizontal flows are stored on the edges of the UgridDataArray topology. The other flows are generally stored on the faces; this includes the flow-lower-face. We'll create a dataset for the horizontal flows for further analysis. .. GENERATED FROM PYTHON SOURCE LINES 201-207 .. code-block:: Python cbc_grid = cbc["flow-horizontal-face-x"].grid ds = xu.UgridDataset(grids=cbc_grid) ds["u"] = cbc["flow-horizontal-face-x"] ds["v"] = cbc["flow-horizontal-face-y"] .. GENERATED FROM PYTHON SOURCE LINES 208-214 Visualize the results --------------------- We can quickly and easily visualize the output with the plotting functions provided by xarray and xugrid. We'll add some some edge coordinates to the dataset so that they can be used to place the arrows in the quiver plot. .. GENERATED FROM PYTHON SOURCE LINES 214-223 .. code-block:: Python ds = ds.ugrid.assign_edge_coords() fig, ax = plt.subplots() head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax) ds.isel(time=0, layer=0).plot.quiver( x="mesh2d_edge_x", y="mesh2d_edge_y", u="u", v="v", color="white" ) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_003.png :alt: layer = 1, time = 1.0 :srcset: /examples/mf6/images/sphx_glr_circle_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 224-226 As would be expected from our model input, we observe circular groundwater mounding and increasing flows as we move from the center to the exterior. .. GENERATED FROM PYTHON SOURCE LINES 228-235 Slice the model domain ---------------------- We may also quickly setup a smaller model. We'll select half of the original domain. To set up the boundary conditions on the clipped edges you can provide a states_for_boundary dictionary. In this case we add the head values of the computed full domain simulation as the clipped boundary values .. GENERATED FROM PYTHON SOURCE LINES 235-244 .. code-block:: Python states_for_boundary = { "GWF_1": head.compute(), } half_simulation = simulation.clip_box( x_max=0.0, states_for_boundary=states_for_boundary ) .. GENERATED FROM PYTHON SOURCE LINES 245-246 Let's run the model, read the results, and visualize. .. GENERATED FROM PYTHON SOURCE LINES 246-252 .. code-block:: Python modeldir = imod.util.temporary_directory() half_simulation.write(modeldir) half_simulation.run() head = half_simulation.open_head() .. GENERATED FROM PYTHON SOURCE LINES 253-254 Let's add constant head boundaries together and plot them .. GENERATED FROM PYTHON SOURCE LINES 254-270 .. code-block:: Python half_simulation_constant_head = half_simulation["GWF_1"]["chd"]["head"] clipped_half_simulation_constant_head = ( half_simulation["GWF_1"]["chd_clipped"]["head"].sel(layer=2).isel(time=0) ) all_boundaries_constant_head = half_simulation_constant_head.where( ~isnull(half_simulation_constant_head), clipped_half_simulation_constant_head ) # plot boundary conditions fig, ax = plt.subplots() all_boundaries_constant_head.ugrid.plot(ax=ax) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_004.png :alt: layer = 2, time = 1.0 :srcset: /examples/mf6/images/sphx_glr_circle_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 271-272 plot computed heads .. GENERATED FROM PYTHON SOURCE LINES 272-275 .. code-block:: Python fig, ax = plt.subplots() head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_005.png :alt: layer = 1, time = 1.0 :srcset: /examples/mf6/images/sphx_glr_circle_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.243 seconds) .. _sphx_glr_download_examples_mf6_circle.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: circle.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: circle.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: circle.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_