.. 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:: default .. GENERATED FROM PYTHON SOURCE LINES 22-23 We'll start with the following imports: .. GENERATED FROM PYTHON SOURCE LINES 23-31 .. code-block:: default import matplotlib.pyplot as plt import numpy as np import xarray as xr import xugrid as xu import imod .. GENERATED FROM PYTHON SOURCE LINES 32-40 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 40-46 .. code-block:: default grid = imod.data.circle() grid .. rst-class:: sphx-glr-script-out .. code-block:: none Dimensions: (mesh2d_nFaces: 216, mesh2d_nMax_face_nodes: 3, mesh2d_nEdges: 342, two: 2, mesh2d_nNodes: 127) Coordinates: mesh2d_node_x (mesh2d_nNodes) float64 0.0 6.123e-14 ... 331.1 157.5 mesh2d_node_y (mesh2d_nNodes) float64 0.0 1e+03 500.0 ... 764.7 818.3 Dimensions without coordinates: mesh2d_nFaces, mesh2d_nMax_face_nodes, mesh2d_nEdges, two, mesh2d_nNodes Data variables: mesh2d int64 0 mesh2d_face_nodes (mesh2d_nFaces, mesh2d_nMax_face_nodes) int32 0 7 ... 1 mesh2d_edge_nodes (mesh2d_nEdges, two) int64 0 7 0 17 0 ... 124 126 125 126 Attributes: Conventions: CF-1.8 UGRID-1.0 .. GENERATED FROM PYTHON SOURCE LINES 47-48 We can plot this object as follows: .. GENERATED FROM PYTHON SOURCE LINES 48-53 .. code-block:: default 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 54-70 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 70-88 .. code-block:: default 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 89-92 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 92-103 .. code-block:: default 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 104-109 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 109-153 .. code-block:: default 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(head=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", csv_output=False, no_ptc=True, 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 154-155 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 155-159 .. code-block:: default modeldir = imod.util.temporary_directory() simulation.write(modeldir) .. rst-class:: sphx-glr-script-out .. code-block:: none /builds/deltares/imod/imod-python/imod/mf6/disv.py:126: FutureWarning: the 'line_terminator'' keyword is deprecated, use 'lineterminator' instead. self._verts_dataframe().to_csv( /builds/deltares/imod/imod-python/imod/mf6/disv.py:132: FutureWarning: the 'line_terminator'' keyword is deprecated, use 'lineterminator' instead. self._cell2d_dataframe().to_csv( .. GENERATED FROM PYTHON SOURCE LINES 160-168 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 168-171 .. code-block:: default simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 172-176 Open the results ---------------- First, we'll open the heads (.hds) file. .. GENERATED FROM PYTHON SOURCE LINES 176-183 .. code-block:: default head = imod.mf6.out.open_hds( modeldir / "GWF_1/GWF_1.hds", modeldir / "GWF_1/disv.disv.grb", ) head .. raw:: html 
<xarray.DataArray 'head' (time: 1, layer: 2, mesh2d_nFaces: 216)>
    dask.array<stack, shape=(1, 2, 216), dtype=float64, chunksize=(1, 2, 216), chunktype=numpy.ndarray>
    Coordinates:
      * layer          (layer) int64 1 2
      * time           (time) float64 1.0
      * mesh2d_nFaces  (mesh2d_nFaces) int64 0 1 2 3 4 5 ... 210 211 212 213 214 215


.. GENERATED FROM PYTHON SOURCE LINES 184-189 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 189-196 .. code-block:: default cbc = imod.mf6.open_cbc( modeldir / "GWF_1/GWF_1.cbc", modeldir / "GWF_1/disv.disv.grb", ) print(cbc.keys()) .. rst-class:: sphx-glr-script-out .. code-block:: none dict_keys(['flow-horizontal-face', 'flow-horizontal-face-x', 'flow-horizontal-face-y', 'flow-lower-face', 'chd']) .. GENERATED FROM PYTHON SOURCE LINES 197-206 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 206-211 .. code-block:: default ds = xu.UgridDataset(grids=grid) ds["u"] = cbc["flow-horizontal-face-x"] ds["v"] = cbc["flow-horizontal-face-y"] .. GENERATED FROM PYTHON SOURCE LINES 212-218 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 218-227 .. code-block:: default 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 228-230 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. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.763 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-python :download:`Download Python source code: circle.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: circle.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_