quickstart is so close

Author: wpbonelli

docs/examples/quickstart.py

@@ -5,23 +5,32 @@
 from flopy.discretization.modeltime import ModelTime
 from flopy.discretization.structuredgrid import StructuredGrid
 
-from flopy4.mf6.gwf import Chd, Gwf, Npf, Oc
+from flopy4.mf6.gwf import Chd, Gwf, Ic, Npf, Oc
 from flopy4.mf6.ims import Ims
 from flopy4.mf6.simulation import Simulation
 
 name = "quickstart"
 workspace = Path(__file__).parent / name
 time = ModelTime(perlen=[1.0], nstp=[1])
-grid = StructuredGrid(nlay=1, nrow=10, ncol=10)
+grid = StructuredGrid(
+    nlay=1,
+    nrow=10,
+    ncol=10,
+    delr=1.0 * np.ones(10),
+    delc=1.0 * np.ones(10),
+    top=1.0 * np.ones((10, 10)),
+    botm=0.0 * np.ones((1, 10, 10)),
+)
 sim = Simulation(name=name, workspace=workspace, tdis=time)
-ims = Ims(parent=sim)
 gwf_name = "mymodel"
+ims = Ims(parent=sim, models=[gwf_name])  # temporary hack
 gwf = Gwf(parent=sim, name=gwf_name, save_flows=True, dis=grid)
 npf = Npf(parent=gwf, save_specific_discharge=True)
 chd = Chd(
     parent=gwf,
     head={0: {(0, 0, 0): 1.0, (0, 9, 9): 0.0}},
 )
+ic = Ic(parent=gwf, strt=1.0)
 oc = Oc(
     parent=gwf,
     budget_file=f"{gwf.name}.bud",
@@ -30,18 +39,15 @@
     save_budget={0: "all"},
 )
 
-# sim.write()
+sim.write()
 sim.run(verbose=True)
 
-# check CHD
 assert chd.data["head"][0, 0] == 1.0
 assert chd.data.head.sel(per=0)[99] == 0.0
 assert np.allclose(chd.data.head[:, 1:99], np.full(98, 1e30))
 
-# check DIS
 assert gwf.dis.data.botm.sel(lay=0, col=0, row=0) == 0.0
 
-# check OC
 assert oc.data["save_head"][0] == "all"
 assert oc.data.save_head.sel(per=0) == "all"
 

flopy4/mf6/converter.py

@@ -13,6 +13,7 @@
 from flopy4.mf6.exchange import Exchange
 from flopy4.mf6.model import Model
 from flopy4.mf6.package import Package
+from flopy4.mf6.solution import Solution
 from flopy4.mf6.spec import fields_dict, get_blocks
 
 
@@ -45,9 +46,10 @@ def _get_binding_type(component: Component) -> str:
         def _get_binding_terms(component: Component) -> tuple[str, ...] | None:
             if isinstance(component, Exchange):
                 return (component.exgmnamea, component.exgmnameb)  # type: ignore
+            elif isinstance(component, Solution):
+                return tuple(component.models)
             elif isinstance(component, (Model, Package)):
                 return (component.name,)  # type: ignore
-            # TODO solutions
             return None
 
         return cls(
@@ -178,7 +180,18 @@ def unstructure_component(value: Component) -> dict[str, Any]:
             _attach_field_metadata(dataset, type(value), list(block.keys()))
             blocks[f"{block_name} {kper + 1}"] = {block_name: dataset}
 
-    return {name: block for name, block in blocks.items() if block}
+    # make sure options block always comes first
+    if "options" in blocks:
+        options_block = blocks.pop("options")
+        blocks = {"options": options_block, **blocks}
+
+    # total temporary hack! manually set solutiongroup 1. still need to support multiple..
+    if "solutiongroup" in blocks:
+        sg = blocks["solutiongroup"]
+        blocks["solutiongroup 1"] = sg
+        del blocks["solutiongroup"]
+
+    return {name: block for name, block in blocks.items() if name != "period"}
 
 
 def _make_converter() -> Converter:

flopy4/mf6/gwf/__init__.py

@@ -31,6 +31,11 @@ def convert_grid(value):
 
 @xattree
 class Gwf(Model):
+    @define
+    class NewtonOptions:
+        newton: bool = field()
+        under_relaxation: bool = field()
+
     @define
     class Output:
         parent: "Gwf" = attrs.field(repr=False)
@@ -51,6 +56,14 @@ def budget(self):
                 self.parent.parent.workspace / f"{self.parent.name}.dis.grb",
             )
 
+    _list: Optional[str] = field(block="options", default=None)
+    print_input: bool = field(block="options", default=False)
+    print_flows: bool = field(block="options", default=False)
+    save_flows: bool = field(block="options", default=False)
+    newtonoptions: Optional[NewtonOptions] = field(block="options", default=None)
+    nc_mesh2d_filerecord: Optional[Path] = field(block="options", default=None)
+    nc_structured_filerecord: Optional[Path] = field(block="options", default=None)
+    nc_filerecord: Optional[Path] = field(block="options", default=None)
     dis: Dis = field(converter=convert_grid, block="packages")
     ic: Ic = field(block="packages")
     oc: Oc = field(block="packages")
@@ -62,20 +75,6 @@ def budget(self):
         default=attrs.Factory(lambda self: Gwf.Output(self), takes_self=True)
     )
 
-    @define
-    class NewtonOptions:
-        newton: bool = field()
-        under_relaxation: bool = field()
-
-    list: Optional[str] = field(block="options", default=None)
-    print_input: bool = field(block="options", default=False)
-    print_flows: bool = field(block="options", default=False)
-    save_flows: bool = field(block="options", default=False)
-    newtonoptions: Optional[NewtonOptions] = field(block="options", default=None)
-    nc_mesh2d_filerecord: Optional[Path] = field(block="options", default=None)
-    nc_structured_filerecord: Optional[Path] = field(block="options", default=None)
-    nc_filerecord: Optional[Path] = field(block="options", default=None)
-
     @property
     def grid(self) -> Grid:
         return self.dis.to_grid()

flopy4/mf6/gwf/ic.py

@@ -8,13 +8,13 @@
 from flopy4.mf6.spec import array, field
 
 
-@xattree
+@xattree(kw_only=True)
 class Ic(Package):
+    export_array_ascii: bool = field(block="options", default=False)
+    export_array_netcdf: bool = field(block="options", default=False)
     strt: NDArray[np.float64] = array(
-        block="packagedata",
+        block="griddata",
         dims=("nnodes",),
         default=1.0,
         converter=Converter(dict_to_array, takes_self=True, takes_field=True),
     )
-    export_array_ascii: bool = field(block="options", default=False)
-    export_array_netcdf: bool = field(block="options", default=False)

flopy4/mf6/ims.py

@@ -12,7 +12,7 @@
 class Ims(Solution):
     solution_package: ClassVar[Sln] = Sln(abbr="ims", pattern="*")
 
-    print_option: bool = field(block="options", default=False)
+    print_option: Optional[str] = field(block="options", default=None)
     complexity: str = field(block="options", default="simple")
     csv_outer_output_file: Optional[Path] = field(default=None, block="options")
     csv_inner_output_file: Optional[Path] = field(block="options", default=None)

flopy4/mf6/simulation.py

@@ -23,10 +23,10 @@ def convert_time(value):
 
 @xattree
 class Simulation(Context):
+    tdis: Tdis = field(converter=convert_time, block="timing")
     models: dict[str, Model] = field(block="models")
     exchanges: dict[str, Exchange] = field(block="exchanges")
     solutions: dict[str, Solution] = field(block="solutiongroup")
-    tdis: Tdis = field(converter=convert_time, block="timing")
     filename: str = field(default="mfsim.nam", init=False)
 
     def __attrs_post_init__(self):

flopy4/mf6/solution.py

@@ -1,10 +1,11 @@
 from abc import ABC
 
+import attrs
 from xattree import xattree
 
 from flopy4.mf6.package import Package
 
 
 @xattree
 class Solution(Package, ABC):
-    pass
+    models: list[str] = attrs.field()