baseline twri example

Author: mjreno

docs/examples/image/qs.png

@@ -11,6 +11,7 @@
 
 name = "quickstart"
 workspace = Path(__file__).parent / name
+workspace.mkdir(exist_ok=True)
 time = Time(perlen=[1.0], nstp=[1])
 grid = StructuredGrid(
     nlay=1,
@@ -65,5 +66,5 @@
 ax.grid(which="both", color="white")
 head.plot.imshow(ax=ax)
 head.plot.contour(ax=ax, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
-budget.plot.quiver(x="x", y="y", u="npf-qx", v="npf-qy", ax=ax, color="white")
+# budget.plot.quiver(x="x", y="y", u="npf-qx", v="npf-qy", ax=ax, color="white")
 fig.savefig(workspace / ".." / "image" / "quickstart.png")

docs/examples/image/qsprj.png

@@ -5,7 +5,6 @@
 from pathlib import Path
 
 import numpy as np
-import xarray as xr
 
 import flopy4
 
@@ -16,7 +15,7 @@
 nper = time.nper
 
 # Grid
-nlay = 5
+nlay = 3
 nrow = 15
 ncol = 15
 shape = (nlay, nrow, ncol)
@@ -25,18 +24,19 @@
 delc = np.ones((nrow,)) * 5000.0
 idomain = np.ones(shape, dtype=int)
 top = np.ones((nrow, ncol), dtype=float) * 200.0
-bottom = np.stack([np.full((nrow, ncol), val) for val in [-150.0, -200.0, -300.0, -350.0, -450.0]])
+bottom = np.stack([np.full((nrow, ncol), val) for val in [-200.0, -300.0, -450.0]])
 grid = flopy4.mf6.utils.grid.StructuredGrid(
     nlay=nlay, nrow=nrow, ncol=ncol, top=top, botm=bottom, delr=delr, delc=delc, idomain=idomain
 )
 dims = {"nper": nper, **dict(grid.dataset.sizes)}  # TODO: temporary
 
 # Grid discretization
+# TODO: xorigin, yorigin
 dis = flopy4.mf6.gwf.Dis.from_grid(grid=grid)
 
 # Constant head boundary on the left
 chd = flopy4.mf6.gwf.Chd(
-    head={"*": {(k, i, 0): 0.0 for k in range(nlay) for i in range(nrow)}},
+    head={"*": {(k, i, 0): 0.0 for k in range(nlay - 1) for i in range(nrow)}},
     print_input=True,
     print_flows=True,
     save_flows=True,
@@ -47,8 +47,8 @@
 elevation = [0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0]
 conductance = 1.0
 drn = flopy4.mf6.gwf.Drn(
-    elev={"*": {(0, 7, j): elevation[j] for j in range(9)}},
-    cond={"*": {(0, 7, j): conductance for j in range(9)}},
+    elev={"*": {(0, 7, j + 1): elevation[j] for j in range(9)}},
+    cond={"*": {(0, 7, j + 1): conductance for j in range(9)}},
     print_input=True,
     print_flows=True,
     save_flows=True,
@@ -59,15 +59,15 @@
 ic = flopy4.mf6.gwf.Ic(strt=0.0, dims=dims)
 
 # Node properties
-icelltype = np.stack([np.full((nrow, ncol), val) for val in [1, 0, 0, 0, 0]])
-k = np.stack([np.full((nrow, ncol), val) for val in [1.0e-3, 1.0e-8, 1.0e-4, 5.0e-7, 2.0e-4]])
-k33 = np.stack([np.full((nrow, ncol), val) for val in [1.0e-3, 1.0e-8, 1.0e-4, 5.0e-7, 2.0e-4]])
+icelltype = np.stack([np.full((nrow, ncol), val) for val in [1, 0, 0]])
+k = np.stack([np.full((nrow, ncol), val) for val in [1.0e-3, 1.0e-4, 2.0e-4]])
+k33 = np.stack([np.full((nrow, ncol), val) for val in [2.0e-8, 2.0e-8, 2.0e-8]])
 npf = flopy4.mf6.gwf.Npf(
     # TODO: no need for reshaping once array structuring converter is done
     icelltype=icelltype.reshape((nodes,)),
     k=k.reshape((nodes,)),
     k33=k33.reshape((nodes,)),
-    cvoptions=flopy4.mf6.gwf.Npf.CvOptions(variablecv=True, dewatered=True),
+    cvoptions=flopy4.mf6.gwf.Npf.CvOptions(dewatered=True),
     perched=True,
     save_flows=True,
     dims=dims,
@@ -78,26 +78,34 @@
     storagecoefficient=False,
     ss=1.0e-5,
     sy=0.15,
-    transient={"*": False},
+    steady_state=[True, False, False],
     iconvert=0,
     dims=dims,
 )
 
 # Uniform recharge on the top layer
-rch_rate = xr.full_like(grid.idomain.sel(k=1), 3.0e-8, dtype=float)
-rch = flopy4.mf6.gwf.Rch(recharge=rch_rate, dims=dims)
+rch_rate = np.stack(np.full((nlay, nrow, ncol), flopy4.mf6.constants.FILL_DNODATA))
+rate = np.repeat(np.expand_dims(rch_rate, axis=0), repeats=nper, axis=0)
+rate[0, 0, :] = 3.0e-8
+rch = flopy4.mf6.gwf.Rch(recharge=rate.reshape(nper, -1), dims=dims)
 
 # Output control
 # TODO: show both ways to set up the Oc package, strings
 # and proper record types? and/or with perioddata param?
-oc = flopy4.mf6.gwf.Oc(save_head="all", save_budget="all", dims=dims)
+oc = flopy4.mf6.gwf.Oc(
+    budget_file="gwf.bud",
+    head_file="gwf.hds",
+    save_head={0: "all"},
+    save_budget={0: "all"},
+    dims=dims,
+)
 
 # Wells scattered throughout the model
 wel_q = -5.0
 wel_nodes = [
-    [3, 4, 10, -5.0],
-    [2, 3, 5, -5.0],
-    [2, 5, 11, -5.0],
+    [2, 4, 10, -5.0],
+    [1, 3, 5, -5.0],
+    [1, 5, 11, -5.0],
     [0, 8, 7, -5.0],
     [0, 8, 9, -5.0],
     [0, 8, 11, -5.0],
@@ -119,14 +127,14 @@
 # Flow model
 gwf = flopy4.mf6.gwf.Gwf(
     dis=grid,
-    chd=chd,
     ic=ic,
     npf=npf,
     sto=sto,
     oc=oc,
-    rch=rch,
-    wel=wel,
-    drn=drn,
+    chd=[chd],
+    rch=[rch],
+    wel=[wel],
+    drn=[drn],
     dims=dims,
 )
 
@@ -143,6 +151,7 @@
     scaling_method=None,
     reordering_method=None,
     relaxation_factor=0.97,
+    models=["gwf"],
 )
 
 # TDIS
@@ -166,10 +175,9 @@
 sim.run()  # assumes the ``mf6`` executable is available on your PATH.
 
 # Load head results
-gwf_ws = Path(workspace) / gwf.name
 head = flopy4.mf6.utils.open_hds(
-    gwf_ws / f"{gwf.name}.hds",
-    gwf_ws / f"{dis.name}.dis.grb",
+    workspace / f"{gwf.name}.hds",
+    workspace / f"{gwf.name}.dis.grb",
 )
 
 # Plot head results

docs/examples/image/quickstart.png

@@ -107,7 +107,7 @@ def array2string(value: NDArray) -> str:
     format = (
         "%d"
         if np.issubdtype(value.dtype, np.integer)
-        else "%f"
+        else "%.9f"
         if np.issubdtype(value.dtype, np.floating)
         else "%s"
     )

docs/examples/qsplot.py

@@ -112,12 +112,18 @@ def oc_setting_data(rec):
     match value.dtype:
         case np.bool:
             # supports boolean dataarrays, e.g. STO steady_state and transient
-            dat = {kper: "" for kper in range(value.sizes["nper"]) if value.values[kper]}  # type: ignore
-            data[name] = dat
+            # e.g. steady_state to steady-state, why is't this the dfn name?
+            fname = name.replace("_", "-")  # type: ignore
+            dat = {kper: "" for kper in range(value.sizes["nper"]) if value.values[kper]}
+            data[fname] = dat
         case np.dtypes.StringDType():
             # supports string dataarrays, e.g. OC save_budget, save_head
             fname = name.replace("_", " ")
-            dat = {kper: value.values[kper] for kper in range(value.sizes["nper"])}
+            dat = {
+                kper: value.values[kper]
+                for kper in range(value.sizes["nper"])
+                if value.values[kper].lower() in ["first", "last", "steps", "all"]
+            }
             data[fname] = dat
         case object():
             # supports object dataararys, e.g. OC PrintSaveSetting
@@ -139,6 +145,8 @@ def oc_setting_data(rec):
 
 
 def unstructure_component(value: Component) -> dict[str, Any]:
+    from flopy4.mf6.constants import FILL_DNODATA
+
     blockspec = dict(sorted(value.dfn.blocks.items(), key=block_sort_key))  # type: ignore
     blocks: dict[str, dict[str, Any]] = {}
     xatspec = xattree.get_xatspec(type(value))
@@ -236,18 +244,21 @@ def unstructure_component(value: Component) -> dict[str, Any]:
 
         # setup indexed period blocks, combine arrays into datasets
         for kper, block in period_blocks.items():
-            blocks[f"period {kper + 1}"] = {}
+            key = f"period {kper + 1}"
             for arr_name, val in block.items():
-                match block[arr_name]:
-                    case str():
-                        # non data period parameters have their period
-                        # write key set in the _hack_period_non_numeric
-                        # routine
-                        blocks[f"period {kper + 1}"][arr_name] = val
-                    case xr.DataArray():
-                        blocks[f"period {kper + 1}"]["period"] = xr.Dataset(
-                            block, coords=block[arr_name].coords
-                        )
+                if np.any(val != FILL_DNODATA):
+                    if key not in blocks:
+                        blocks[key] = {}
+                    match block[arr_name]:
+                        case str():
+                            # non data period parameters have their period
+                            # write key set in the _hack_period_non_numeric
+                            # routine
+                            blocks[f"period {kper + 1}"][arr_name] = val
+                        case xr.DataArray():
+                            blocks[f"period {kper + 1}"]["period"] = xr.Dataset(
+                                block, coords=block[arr_name].coords
+                            )
 
         # combine "perioddata" block arrays (tdis, ats) into datasets
         # so they render as lists. temp hack TODO do this generically

docs/examples/quickstart.py

@@ -73,7 +73,7 @@ def test_dumps_sto():
     print("STO dump:")
     print(dumped)
     assert "BEGIN PERIOD 1\n TRANSIENT" in dumped
-    assert "BEGIN PERIOD 2\n STEADY_STATE" in dumped
+    assert "BEGIN PERIOD 2\n STEADY-STATE" in dumped
     assert "BEGIN PERIOD 3\n TRANSIENT" in dumped
     assert dumped