reshaping no longer necessary (#279)

get rid of reshaping in TWRI even for grid and layer array package flavors

Author: wpbonelli

docs/dev/array-converter-design.md

@@ -251,10 +251,15 @@ def _reshape_grid(
     """
     Perform structured↔flat grid conversion.
 
-    Handles:
+    Handles full 3D grids (nodes dimension):
     - (nlay, nrow, ncol) → (nodes,)
     - (nper, nlay, nrow, ncol) → (nper, nodes)
-    - Preserves xarray metadata if input is xarray
+
+    Handles per-layer 2D arrays (ncpl dimension):
+    - (nrow, ncol) → (ncpl,)
+    - (nper, nrow, ncol) → (nper, ncpl)
+
+    Preserves xarray metadata if input is xarray
     """
 ```
 
@@ -539,19 +544,54 @@ coords = {coord_name: self._coords[coord_name], "x": self._coords["x"]}
 
 **Location**: `flopy4/mf6/utils/grid.py` lines 261-284
 
+#### 4. Missing Converters on Grid-Based Packages
+**Issue**: Grid-based ("g") and array-based ("a") package variants (Chdg, Welg, Drng, Rcha) were missing `converter=Converter(structure_array, ...)` on their period block array fields. Without converters, xattree validates dimension counts before reshape logic runs, causing errors when passing structured arrays.
+
+**Solution**: Added converters to all grid-based and array-based package fields:
+- `Chdg.head`: Added converter for automatic reshaping
+- `Welg.q`: Added converter for automatic reshaping
+- `Drng.elev` and `Drng.cond`: Added converters for automatic reshaping
+- `Rcha.recharge`: Added converter for automatic reshaping
+
+**Location**: `flopy4/mf6/gwf/chdg.py`, `welg.py`, `drng.py`, `rcha.py`
+
+#### 5. Per-Layer Array Dimension (ncpl)
+**Issue**: Array-based packages use `ncpl` ("number of cells per layer") dimension for 2D per-layer data. The reshape detection only handled `nodes` (3D full grids) but not `ncpl` (2D per-layer), causing errors like:
+```
+ValueError: Shape mismatch: (3, 15, 15) vs (3, 225)
+```
+
+**Solution**: Extended `_detect_grid_reshape` to handle `ncpl` dimension:
+```python
+# Handle 'ncpl' dimension (cells per layer, 2D per-layer arrays)
+if "ncpl" in expected_dims and has_structured_2d:
+    # Case: (nrow, ncol) → (ncpl,)
+    # Case: (nper, nrow, ncol) → (nper, ncpl)
+```
+
+This allows automatic reshaping for array-based packages like Rcha.
+
+**Location**: `flopy4/mf6/converter/structure.py` lines 113-127
+
 ### Validation Benefits
 
 The stricter validation in the new converter caught the StructuredGrid bug that had existed undetected. While this temporarily broke tests, it exposed a real issue that would have caused problems downstream. This demonstrates the value of proper input validation.
 
 ### User-Facing Improvements
 
 Users can now:
-1. Pass structured arrays `(nlay, nrow, ncol)` directly - automatic reshaping to `(nodes,)`
+1. Pass structured arrays directly - automatic reshaping for all dimensions:
+   - 3D full grids: `(nlay, nrow, ncol)` → `(nodes,)`
+   - 2D per-layer: `(nrow, ncol)` → `(ncpl,)`
+   - Time-varying 3D: `(nper, nlay, nrow, ncol)` → `(nper, nodes)`
+   - Time-varying 2D: `(nper, nrow, ncol)` → `(nper, ncpl)`
 2. Use DataFrames from `package.stress_period_data` to initialize new packages
 3. Mix value types within dicts (scalars, arrays, xarrays, DataFrames)
 4. Rely on strict validation to catch dimension mismatches early
 
-Example - no manual reshaping needed:
+Examples - no manual reshaping needed:
+
+**Example 1: Standard packages (3D grid arrays)**
 ```python
 # Before (manual reshape required)
 icelltype = np.stack([np.full((nrow, ncol), val) for val in [1, 0, 0]])
@@ -561,3 +601,29 @@ npf = Npf(icelltype=icelltype.reshape((nodes,)), ...)
 icelltype = np.stack([np.full((nrow, ncol), val) for val in [1, 0, 0]])
 npf = Npf(icelltype=icelltype, ...)  # Automatically reshaped to (nodes,)
 ```
+
+**Example 2: Grid-based packages (time-varying 3D arrays)**
+```python
+# Before (manual reshape required)
+head = np.full((nper, nlay, nrow, ncol), FILL_DNODATA)
+head[0, :2, :, 0] = 0.0
+chdg = Chdg(head=head.reshape(nper, -1), ...)
+
+# After (automatic reshape)
+head = np.full((nper, nlay, nrow, ncol), FILL_DNODATA)
+head[0, :2, :, 0] = 0.0
+chdg = Chdg(head=head, ...)  # Automatically reshaped to (nper, nodes)
+```
+
+**Example 3: Array-based packages (time-varying 2D per-layer arrays)**
+```python
+# Before (manual reshape required)
+recharge = np.full((nper, nrow, ncol), FILL_DNODATA)
+recharge[0, ...] = 3.0e-8
+rcha = Rcha(recharge=recharge.reshape(nper, -1), ...)
+
+# After (automatic reshape)
+recharge = np.full((nper, nrow, ncol), FILL_DNODATA)
+recharge[0, ...] = 3.0e-8
+rcha = Rcha(recharge=recharge, ...)  # Automatically reshaped to (nper, ncpl)
+```

docs/examples/twri.py

@@ -101,7 +101,7 @@ def plot_head(head, workspace):
 rch_rate = 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)
+rch = flopy4.mf6.gwf.Rch(recharge=rate, dims=dims)
 
 # Output control
 # TODO: show both ways to set up the Oc package, strings
@@ -211,7 +211,7 @@ def plot_head(head, workspace):
     print_input=True,
     print_flows=True,
     save_flows=True,
-    head=head.reshape(nper, -1),
+    head=head,
     dims=dims,
 )
 
@@ -225,8 +225,8 @@ def plot_head(head, workspace):
     print_input=True,
     print_flows=True,
     save_flows=True,
-    elev=elev.reshape(nper, -1),
-    cond=cond.reshape(nper, -1),
+    elev=elev,
+    cond=cond,
     dims=dims,
 )
 
@@ -235,14 +235,14 @@ def plot_head(head, workspace):
 for layer, row, col in wel_nodes:
     q[0, layer, row, col] = wel_q
 welg = flopy4.mf6.gwf.Welg(
-    q=q.reshape(nper, -1),
+    q=q,
     dims=dims,
 )
 
 # recharge
 recharge = np.repeat(np.expand_dims(LAYER_NODATA, axis=0), repeats=nper, axis=0)
 recharge[0, ...] = 3.0e-8
-rcha = flopy4.mf6.gwf.Rcha(recharge=recharge.reshape(nper, -1), dims=dims)
+rcha = flopy4.mf6.gwf.Rcha(recharge=recharge, dims=dims)
 
 # remove list based inputs
 # TODO: show variations on removing packages

flopy4/mf6/converter/structure.py

@@ -86,34 +86,45 @@ def _detect_grid_reshape(
     target_shape : tuple | None
         Target shape for reshape, or None
     """
-    # Check if we expect flat 'nodes' dimension
-    if "nodes" not in expected_dims:
-        return False, None
-
     # Get expected shape
     expected_shape = tuple(dim_dict.get(d, d) for d in expected_dims)
 
     # Check if value has structured dimensions
-    has_structured = "nlay" in dim_dict and "nrow" in dim_dict and "ncol" in dim_dict
-
-    if not has_structured:
-        return False, None
-
-    nlay = dim_dict["nlay"]
-    nrow = dim_dict["nrow"]
-    ncol = dim_dict["ncol"]
-    nodes = dim_dict.get("nodes", nlay * nrow * ncol)
-
-    # Check for structured→flat conversion
-    # Case 1: (nlay, nrow, ncol) → (nodes,)
-    if value_shape == (nlay, nrow, ncol) and expected_shape == (nodes,):
-        return True, (nodes,)
-
-    # Case 2: (nper, nlay, nrow, ncol) → (nper, nodes)
-    if "nper" in expected_dims:
-        nper = dim_dict["nper"]
-        if value_shape == (nper, nlay, nrow, ncol) and expected_shape == (nper, nodes):
-            return True, (nper, nodes)
+    has_structured_3d = "nlay" in dim_dict and "nrow" in dim_dict and "ncol" in dim_dict
+    has_structured_2d = "nrow" in dim_dict and "ncol" in dim_dict
+
+    # Handle 'nodes' dimension (full 3D grid)
+    if "nodes" in expected_dims and has_structured_3d:
+        nlay = dim_dict["nlay"]
+        nrow = dim_dict["nrow"]
+        ncol = dim_dict["ncol"]
+        nodes = dim_dict.get("nodes", nlay * nrow * ncol)
+
+        # Case 1: (nlay, nrow, ncol) → (nodes,)
+        if value_shape == (nlay, nrow, ncol) and expected_shape == (nodes,):
+            return True, (nodes,)
+
+        # Case 2: (nper, nlay, nrow, ncol) → (nper, nodes)
+        if "nper" in expected_dims:
+            nper = dim_dict["nper"]
+            if value_shape == (nper, nlay, nrow, ncol) and expected_shape == (nper, nodes):
+                return True, (nper, nodes)
+
+    # Handle 'ncpl' dimension (cells per layer, 2D per-layer arrays)
+    if "ncpl" in expected_dims and has_structured_2d:
+        nrow = dim_dict["nrow"]
+        ncol = dim_dict["ncol"]
+        ncpl = dim_dict.get("ncpl", nrow * ncol)
+
+        # Case 3: (nrow, ncol) → (ncpl,)
+        if value_shape == (nrow, ncol) and expected_shape == (ncpl,):
+            return True, (ncpl,)
+
+        # Case 4: (nper, nrow, ncol) → (nper, ncpl)
+        if "nper" in expected_dims:
+            nper = dim_dict["nper"]
+            if value_shape == (nper, nrow, ncol) and expected_shape == (nper, ncpl):
+                return True, (nper, ncpl)
 
     return False, None
 
@@ -360,7 +371,7 @@ def _parse_dataframe(
             if has_structured:
                 cellid = (int(row["layer"]), int(row["row"]), int(row["col"]))
             else:
-                cellid = (int(row["node"]),) # type: ignore
+                cellid = (int(row["node"]),)  # type: ignore
 
             # Extract field value
             value = row[field_name]

flopy4/mf6/gwf/chdg.py

@@ -1,10 +1,12 @@
 from pathlib import Path
 from typing import ClassVar, Optional
 
+import attrs
 import numpy as np
 from numpy.typing import NDArray
 from xattree import xattree
 
+from flopy4.mf6.converter import structure_array
 from flopy4.mf6.package import Package
 from flopy4.mf6.spec import array, field, path
 from flopy4.mf6.utils.grid import update_maxbound
@@ -33,6 +35,7 @@ class Chdg(Package):
             "nodes",
         ),
         default=None,
+        converter=attrs.Converter(structure_array, takes_self=True, takes_field=True),
         on_setattr=update_maxbound,
     )
     aux: Optional[NDArray[np.float64]] = array(

flopy4/mf6/gwf/drng.py

@@ -1,10 +1,12 @@
 from pathlib import Path
 from typing import ClassVar, Optional
 
+import attrs
 import numpy as np
 from numpy.typing import NDArray
 from xattree import xattree
 
+from flopy4.mf6.converter import structure_array
 from flopy4.mf6.package import Package
 from flopy4.mf6.spec import array, field, path
 from flopy4.mf6.utils.grid import update_maxbound
@@ -34,6 +36,7 @@ class Drng(Package):
             "nodes",
         ),
         default=None,
+        converter=attrs.Converter(structure_array, takes_self=True, takes_field=True),
         on_setattr=update_maxbound,
     )
     cond: Optional[NDArray[np.float64]] = array(
@@ -43,6 +46,7 @@ class Drng(Package):
             "nodes",
         ),
         default=None,
+        converter=attrs.Converter(structure_array, takes_self=True, takes_field=True),
         on_setattr=update_maxbound,
     )
     aux: Optional[NDArray[np.float64]] = array(

flopy4/mf6/gwf/rcha.py

@@ -1,10 +1,12 @@
 from pathlib import Path
 from typing import ClassVar, Optional
 
+import attrs
 import numpy as np
 from numpy.typing import NDArray
 from xattree import xattree
 
+from flopy4.mf6.converter import structure_array
 from flopy4.mf6.package import Package
 from flopy4.mf6.spec import array, field, path
 from flopy4.utils import to_path
@@ -41,6 +43,7 @@ class Rcha(Package):
             "ncpl",
         ),
         default=None,
+        converter=attrs.Converter(structure_array, takes_self=True, takes_field=True),
     )
     aux: Optional[NDArray[np.float64]] = array(
         block="period",

flopy4/mf6/gwf/welg.py

@@ -1,10 +1,12 @@
 from pathlib import Path
 from typing import ClassVar, Optional
 
+import attrs
 import numpy as np
 from numpy.typing import NDArray
 from xattree import xattree
 
+from flopy4.mf6.converter import structure_array
 from flopy4.mf6.package import Package
 from flopy4.mf6.spec import array, field, path
 from flopy4.mf6.utils.grid import update_maxbound
@@ -39,6 +41,7 @@ class Welg(Package):
             "nodes",
         ),
         default=None,
+        converter=attrs.Converter(structure_array, takes_self=True, takes_field=True),
         on_setattr=update_maxbound,
     )
     aux: Optional[NDArray[np.float64]] = array(