fix(io): extending and testing ascii writes

flopy4/mf6/codec/writer/filters.py

@@ -260,9 +260,14 @@ def data2keystring(value: dict | xr.Dataset):
             return
 
         for field_name in value.data_vars.keys():
+            name = (
+                field_name.replace("_", " ").upper()
+                if np.issubdtype(value.data_vars[field_name].dtype, np.str_)
+                else field_name.upper()
+            )
             field_val = value[field_name]
             if hasattr(field_val, "item"):
                 val = field_val.item()
             else:
                 val = field_val
-            yield (field_name.upper(), val)
+            yield (name, val)

flopy4/mf6/codec/writer/templates/macros.jinja

@@ -1,3 +1,5 @@
+{% set inset = "  " %}
+
 {% macro field(name, value) %}
 {% set format = value|field_format %}
 {% if format in ['keyword', 'integer', 'double precision', 'string'] %}
@@ -15,7 +17,7 @@
 
 {% macro scalar(name, value) %}
 {% set format = value|field_format %}
-{% if value is not none %}{{ name.upper() }}{% if format != 'keyword' %} {{ value }}{% endif %}{% endif %}
+{% if value is not none %}{{ inset ~ name.upper() }}{% if format != 'keyword' %} {{ value }}{% endif %}{% endif %}
 {% endmacro %}
 
 {% macro keystring(name, value) %}
@@ -30,7 +32,7 @@
 {{ field_name.upper() }} {{ field(field_value) }}
 {%- endfor %}
 {% else %}
-{{ value|join(" ") }}
+{{ inset ~ value|join(" ") }}
 {%- endif %}
 {% endmacro %}
 
@@ -55,6 +57,6 @@ OPEN/CLOSE {{ value }}
 
 {% macro list(name, value) %}
 {% for row in (value|data2list) %}
-{{ row|join(" ") }}
+{{ inset ~ row|join(" ") }}
 {% endfor %}
 {% endmacro %}

flopy4/mf6/converter.py

@@ -1,4 +1,4 @@
-from collections.abc import MutableMapping
+from collections.abc import Iterable, MutableMapping
 from datetime import datetime
 from pathlib import Path
 from typing import Any
@@ -47,6 +47,8 @@ def _get_binding_type(component: Component) -> str:
             cls_name = component.__class__.__name__
             if isinstance(component, Exchange):
                 return f"{'-'.join([cls_name[:2], cls_name[3:]]).upper()}6"
+            elif isinstance(component, Solution):
+                return f"{component.slntype}6"
             else:
                 return f"{cls_name.upper()}6"
 
@@ -107,7 +109,7 @@ def unstructure_component(value: Component) -> dict[str, Any]:
                         for comp in field_value.values()
                         if comp is not None
                     ]
-                elif isinstance(field_value, (list, tuple)):
+                elif isinstance(field_value, Iterable):
                     components = [
                         _Binding.from_component(comp).to_tuple()
                         for comp in field_value
@@ -170,16 +172,16 @@ def unstructure_component(value: Component) -> dict[str, Any]:
                                 (
                                     field_value.sizes["nper"],
                                     parent.dims["nlay"],
-                                    parent.dims["ncol"],
                                     parent.dims["nrow"],
+                                    parent.dims["ncol"],
                                 )
                             ),
-                            dims=("nper", "nlay", "ncol", "nrow"),
+                            dims=("nper", "nlay", "nrow", "ncol"),
                             coords={
                                 "nper": field_value.coords["nper"],
                                 "nlay": range(parent.dims["nlay"]),
-                                "ncol": range(parent.dims["ncol"]),
                                 "nrow": range(parent.dims["nrow"]),
+                                "ncol": range(parent.dims["ncol"]),
                             },
                             name=field_value.name,
                         )
@@ -189,12 +191,28 @@ def unstructure_component(value: Component) -> dict[str, Any]:
                         for kper in range(field_value.sizes["nper"])
                     }
                 else:
-                    if block_name not in period_data:
-                        period_data[block_name] = {}
-                    period_data[block_name][field_name] = field_value  # type: ignore
+                    if (
+                        # TODO: refactor
+                        # field_name == "save_budget"
+                        # or field_name == "save_head"
+                        # or field_name == "print_budget"
+                        # or field_name == "print_head"
+                        np.issubdtype(field_value.dtype, np.str_)
+                    ):
+                        period_data[field_name] = {
+                            kper: field_value[kper] for kper in range(field_value.sizes["nper"])
+                        }
+                    else:
+                        if block_name not in period_data:
+                            period_data[block_name] = {}
+                        period_data[block_name][field_name] = field_value  # type: ignore
             else:
                 if field_value is not None:
-                    blocks[block_name][field_name] = field_value
+                    if isinstance(field_value, bool):
+                        if field_value:
+                            blocks[block_name][field_name] = field_value
+                    else:
+                        blocks[block_name][field_name] = field_value
 
         if block_name in period_data and isinstance(period_data[block_name], dict):
             dataset = xr.Dataset(period_data[block_name])

flopy4/mf6/ims.py

@@ -11,17 +11,19 @@
 @xattree
 class Ims(Solution):
     solution_package: ClassVar[Sln] = Sln(abbr="ims", pattern="*")
+    slntype: ClassVar[str] = "ims"
 
+    mxiter: Optional[int] = field(default=1)
     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)
     no_ptc: bool = field(default=False, block="options")
     no_ptc_option: Optional[str] = field(default=None, block="options")
     ats_outer_maximum_fraction: Optional[float] = field(block="options", default=None)
-    outer_dvclose: Optional[float] = field(default=None, block="options")
-    outer_maximum: Optional[int] = field(default=None, block="options")
-    under_relaxation: Optional[str] = field(default=None, block="options")
+    outer_dvclose: Optional[float] = field(default=None, block="nonlinear")
+    outer_maximum: Optional[int] = field(default=None, block="nonlinear")
+    under_relaxation: Optional[str] = field(default=None, block="nonlinear")
     under_relaxation_gamma: Optional[float] = field(block="nonlinear", default=None)
     under_relaxation_theta: Optional[float] = field(block="nonlinear", default=None)
     under_relaxation_kappa: Optional[float] = field(block="nonlinear", default=None)

flopy4/mf6/solution.py

@@ -1,11 +1,19 @@
 from abc import ABC
+from pathlib import Path
+from typing import ClassVar, Optional
 
 import attrs
-from xattree import xattree
+from xattree import field, xattree
 
 from flopy4.mf6.package import Package
 
 
 @xattree
 class Solution(Package, ABC):
+    slntype: ClassVar[str] = "sln"
+
+    slnfname: Optional[Path] = field(default=None)  # type: ignore
     models: list[str] = attrs.field(default=attrs.Factory(list))
+
+    def default_filename(self) -> str:
+        return str(self.slnfname) if self.slnfname else f"solution.{self.slntype.lower()}"

test/test_component.py

@@ -274,13 +274,135 @@ def test_ims_dfn():
     assert "inner_maximum" in set(dfn["linear"].keys())
 
 
+def test_gwf_chd01(function_tmpdir):
+    sim_name = "chd01"
+    gwf_name = "gwf_chd01"
+    time = ModelTime(perlen=[5.0], nstp=[1], tsmult=[1.0], time_units="days")
+
+    ims = Ims(
+        slnfname="sln1.ims",
+        models=[gwf_name],
+        print_option="summary",
+        outer_dvclose=1.00000000e-06,
+        outer_maximum=100,
+        under_relaxation="none",
+        inner_maximum=300,
+        inner_dvclose=1.00000000e-06,
+        inner_rclose=1.00000000e-06,
+        linear_acceleration="cg",
+        relaxation_factor=1.0,
+        scaling_method="none",
+        reordering_method="none",
+    )
+
+    sim = Simulation(
+        tdis=time,
+        workspace=function_tmpdir,
+        name=sim_name,
+        solutions={"ims": ims},
+    )
+
+    dis = Dis(
+        nlay=1,
+        nrow=1,
+        ncol=100,
+        delr=1.0,
+        delc=1.0,
+        top=1.0,
+        botm=0.0,
+        idomain=1,
+    )
+
+    gwf = Gwf(parent=sim, save_flows=True, dis=dis, name=gwf_name)
+
+    ic = Ic(parent=gwf, strt=1.0)
+
+    oc = Oc(
+        parent=gwf,
+        budget_file=f"{gwf_name}.cbc",
+        head_file=f"{gwf_name}.hds",
+        # COLUMNS  10  WIDTH  15  DIGITS  6  GENERAL
+        save_head=["last"],
+        # save_head={0: "last"},
+        save_budget=["last"],
+        print_head=["last"],
+        print_budget=["last"],
+    )
+
+    npf = Npf(
+        parent=gwf,
+        save_specific_discharge=True,
+        k=1.0,
+        k33=1.0,
+        icelltype=0,
+    )
+
+    chd = Chd(
+        parent=gwf,
+        print_flows=True,
+        head={0: {(0, 0, 0): 1.0, (0, 0, 99): 0.0}},
+        name="chd-1",
+    )
+
+    sim.write()
+    sim.run()
+
+
+def test_quickstart(function_tmpdir):
+    sim_name = "quickstart"
+    gwf_name = "mymodel"
+    time = ModelTime(perlen=[1.0], nstp=[1], tsmult=[1.0])
+    ims = Ims(models=[gwf_name])
+    dis = Dis(
+        nlay=1,
+        nrow=10,
+        ncol=10,
+        top=1.0,
+        botm=0.0,
+    )
+    sim = Simulation(
+        tdis=time,
+        workspace=function_tmpdir,
+        name=sim_name,
+        solutions={"ims": ims},
+    )
+    gwf = Gwf(parent=sim, dis=dis, name=gwf_name)
+    ic = Ic(parent=gwf)
+    oc = Oc(
+        parent=gwf,
+        budget_file=f"{gwf_name}.bud",
+        head_file=f"{gwf_name}.hds",
+        save_head=["all"],
+        save_budget=["all"],
+    )
+    npf = Npf(parent=gwf, icelltype=0, k=1.0)
+    chd = Chd(parent=gwf, head={0: {(0, 0, 0): 1.0, (0, 9, 9): 0.0}})
+
+    sim.write()
+    sim.run()
+
+
 def test_write_ascii(function_tmpdir):
     sim_name = "sim"
-    time = ModelTime(perlen=[1.0], nstp=[1], tsmult=[1.0])
-    grid = StructuredGrid(nlay=1, nrow=10, ncol=10)
-    sim = Simulation(tdis=time, workspace=function_tmpdir, name=sim_name)
     gwf_name = "gwf"
-    gwf = Gwf(parent=sim, dis=grid, name=gwf_name)
+    time = ModelTime(perlen=[1.0], nstp=[1], tsmult=[1.0])
+    ims = Ims(models=[gwf_name])
+    dis = Dis(
+        nlay=1,
+        nrow=10,
+        ncol=10,
+        delr=1.0,
+        delc=1.0,
+        top=1.0,
+        botm=0.0,
+    )
+    sim = Simulation(
+        tdis=time,
+        workspace=function_tmpdir,
+        name=sim_name,
+        solutions={"ims": ims},
+    )
+    gwf = Gwf(parent=sim, dis=dis, name=gwf_name)
     ic = Ic(parent=gwf)
     oc = Oc(parent=gwf)
     npf = Npf(parent=gwf)

test/test_interface.py

@@ -242,9 +242,6 @@ def test_flopy3_package(tmp_path):
 
 
 def norun_test_flopy3_cbd_small(tmp_path):
-    import sys
-
-    sys.path.append("/home/mjreno/.clone/usgs/flopy/autotest")
     from test_grid_cases import GridCases
 
     time = ModelTime(perlen=[1.0], nstp=[1], tsmult=[1.0])