write runnable gwf base simulation

Author: mjreno

docs/examples/array_example.py

@@ -45,11 +45,12 @@
 constant = data_path / "constant.txt"
 external = data_path / "external.txt"
 shape = (1000, 100)
+type = "double"
 
 # Open and load a NumPy array representation
 
 fhandle = open(internal)
-imfa = MFArray.load(fhandle, data_path, shape, header=False)
+imfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 
 # Get values
 
@@ -70,7 +71,7 @@
 plt.colorbar()
 
 fhandle = open(constant)
-cmfa = MFArray.load(fhandle, data_path, shape, header=False)
+cmfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 cvals = cmfa.value
 plt.imshow(cvals[0:100])
 plt.colorbar()
@@ -93,7 +94,7 @@
 # External
 
 fhandle = open(external)
-emfa = MFArray.load(fhandle, data_path, shape, header=False)
+emfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 evals = emfa.value
 evals
 
@@ -118,7 +119,9 @@
 
 fhandle = open(ilayered)
 shape = (3, 1000, 100)
-ilmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+ilmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 vals = ilmfa.value
 
 ilmfa._value  # internal storage
@@ -165,7 +168,9 @@
 
 fhandle = open(clayered)
 shape = (3, 1000, 100)
-clmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+clmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 
 clmfa._value
 
@@ -218,7 +223,9 @@
 
 fhandle = open(mlayered)
 shape = (3, 1000, 100)
-mlmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+mlmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 
 mlmfa.how
 

flopy4/array.py

@@ -434,6 +434,16 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
         model_shape = kwargs.pop("model_shape", None)
         params = kwargs.pop("blk_params", {})
         mempath = kwargs.pop("mempath", None)
+        atype = kwargs.get("type", None)
+
+        if atype is not None:
+            if atype == "integer":
+                dtype = np.int32
+            elif atype == "double":
+                dtype = np.float64
+        else:
+            raise ValueError("array spec type not defined")
+
         if model_shape and isinstance(shape, str):
             if shape == "(nodes)":
                 n = math.prod([x for x in model_shape])
@@ -459,7 +469,7 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
             lshp = shape[1:]
             objs = []
             for _ in range(nlay):
-                mfa = cls._load(f, cwd, lshp, name)
+                mfa = cls._load(f, cwd, lshp, dtype=dtype, name=name)
                 objs.append(mfa)
 
             return MFArray(
@@ -474,11 +484,17 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
         else:
             kwargs.pop("layered", None)
             return cls._load(
-                f, cwd, shape, layered=layered, name=name, **kwargs
+                f,
+                cwd,
+                shape,
+                layered=layered,
+                dtype=dtype,
+                name=name,
+                **kwargs,
             )
 
     @classmethod
-    def _load(cls, f, cwd, shape, layered=False, **kwargs):
+    def _load(cls, f, cwd, shape, layered=False, dtype=None, **kwargs):
         control_line = multi_line_strip(f).split()
 
         if CommonNames.iprn.lower() in control_line:
@@ -491,25 +507,31 @@ def _load(cls, f, cwd, shape, layered=False, **kwargs):
         clpos = 1
 
         if how == MFArrayType.internal:
-            array = cls.read_array(f)
+            array = cls.read_array(f, dtype)
 
         elif how == MFArrayType.constant:
-            array = float(control_line[clpos])
+            if dtype == np.float64:
+                array = float(control_line[clpos])
+            else:
+                array = int(control_line[clpos])
             clpos += 1
 
         elif how == MFArrayType.external:
             extpath = Path(control_line[clpos])
             fpath = cwd / extpath
             with open(fpath) as foo:
-                array = cls.read_array(foo)
+                array = cls.read_array(foo, dtype)
             clpos += 1
 
         else:
             raise NotImplementedError()
 
         factor = None
         if len(control_line) > 2:
-            factor = float(control_line[clpos + 1])
+            if dtype == np.float64:
+                factor = float(control_line[clpos + 1])
+            else:
+                factor = int(control_line[clpos + 1])
 
         return cls(
             shape,
@@ -521,7 +543,7 @@ def _load(cls, f, cwd, shape, layered=False, **kwargs):
         )
 
     @staticmethod
-    def read_array(f):
+    def read_array(f, dtype):
         """
         Read a MODFLOW 6 array from an open file
         into a flat NumPy array representation.
@@ -538,5 +560,5 @@ def read_array(f):
             astr.append(line)
 
         astr = StringIO(" ".join(astr))
-        array = np.genfromtxt(astr).ravel()
+        array = np.genfromtxt(astr, dtype=dtype).ravel()
         return array

flopy4/block.py

@@ -225,6 +225,7 @@ def load(cls, f, **kwargs):
         name = None
         index = None
         found = False
+        period = False
         params = dict()
         members = cls.params
 
@@ -238,12 +239,17 @@ def load(cls, f, **kwargs):
             if line == "\n":
                 continue
             words = strip(line).lower().split()
-            key = words[0]
+            if period:
+                key = "stress_period_data"
+            else:
+                key = words[0]
             if key == "begin":
                 found = True
                 name = words[1]
                 if len(words) > 2 and str.isdigit(words[2]):
                     index = int(words[2])
+                if name == "period":
+                    period = True
             elif key == "end":
                 break
             elif found:
@@ -268,20 +274,22 @@ def load(cls, f, **kwargs):
                     # TODO: inject from model somehow?
                     # and remove special handling here
                     kwrgs["cwd"] = ""
+                    # kwrgs["type"] = param.type
                     kwrgs["mempath"] = f"{mempath}/{name}"
-                if ptype is not MFArray:
+                if ptype is not MFArray and ptype is not MFList:
                     kwrgs.pop("model_shape", None)
                     kwrgs.pop("blk_params", None)
 
                 params[param.name] = ptype.load(f, **kwrgs)
+                period = False
 
         return cls(name=name, index=index, params=params)
 
     def write(self, f):
         """Write the block to file."""
         index = self.index if self.index is not None else ""
         begin = f"BEGIN {self.name.upper()} {index}\n"
-        end = f"END {self.name.upper()}\n"
+        end = f"END {self.name.upper()}\n\n"
 
         f.write(begin)
         super().write(f)

flopy4/compound.py

@@ -7,7 +7,7 @@
 
 from flopy4.array import MFArray, MFArrayType
 from flopy4.param import MFParam, MFParams, MFReader
-from flopy4.scalar import MFDouble, MFInteger, MFScalar
+from flopy4.scalar import MFScalar
 from flopy4.utils import strip
 
 PAD = "  "
@@ -338,21 +338,39 @@ def load(cls, f, **kwargs) -> "MFList":
         """Load list input with the given component parameters from a file."""
 
         blk_params = kwargs.pop("blk_params", {})
+        model_shape = kwargs.pop("model_shape", None)
         params = kwargs.pop("params", None)
-        type = kwargs.pop("type", None)
         kwargs.pop("mname", None)
-        kwargs.pop("shape", None)
+        kwargs.pop("shape", None)  # e.g. maxbound
+
+        param_lists = []
+        param_cols = []
+        param_types = []
+        for k in list(params):
+            if params[k].name == "aux" or params[k].name == "boundname":
+                continue
+                # raise NotImplementedError(
+                #    "boundames and auxvars not yet supported in period blocks"
+                # )
+            pcols = 0
+            if params[k].shape is None or params[k].shape == "":
+                pcols = 1
+            elif params[k].shape == "(ncelldim)":
+                if model_shape:
+                    pcols = len(model_shape)
+                else:
+                    raise ValueError("model_shape not set")
+            else:
+                pcols = len(params[k].shape.split(","))
+            param_cols.append(pcols)
+            param_lists.append(list())
+            param_types.append(params[k].type)
 
         if list(params.items())[-1][1].shape == "(:)":
-            maxsplit = len(params) - 1
+            maxsplit = sum(param_cols) - 1
         else:
             maxsplit = -1
 
-        param_lists = []
-        # TODO: support multi-dimensional params
-        for i in range(len(params)):
-            param_lists.append(list())
-
         while True:
             pos = f.tell()
             line = f.readline()
@@ -361,9 +379,28 @@ def load(cls, f, **kwargs) -> "MFList":
                 break
             else:
                 tokens = strip(line).split(maxsplit=maxsplit)
-                assert len(tokens) == len(param_lists)
-                for i, token in enumerate(tokens):
-                    param_lists[i].append(token)
+                assert len(tokens) == sum(param_cols)
+                icol = 0
+                for i, p in enumerate(param_lists):
+                    if param_cols[i] == 1:
+                        if param_types[i] == "integer":
+                            param_lists[i].append(int(tokens[icol]))
+                        elif param_types[i] == "double":
+                            param_lists[i].append(float(tokens[icol]))
+                        else:
+                            param_lists[i].append(tokens[icol])
+                        icol += 1
+                    else:
+                        row_l = []
+                        for j in range(param_cols[i]):
+                            if param_types[i] == "integer":
+                                row_l.append(int(tokens[icol]))
+                            elif param_types[i] == "double":
+                                row_l.append(float(tokens[icol]))
+                            else:
+                                row_l.append(tokens[icol])
+                            icol += 1
+                        param_lists[i].append(row_l)
 
         if blk_params and "dimensions" in blk_params:
             nbound = blk_params.get("dimensions").get("nbound")
@@ -372,31 +409,41 @@ def load(cls, f, **kwargs) -> "MFList":
                     if len(param_list) > nbound:
                         raise ValueError("MFList nbound not satisfied")
 
-        list_params = MFList.create_list_params(params, param_lists, **kwargs)
-        return cls(list_params, type=type, **kwargs)
+        list_params = MFList.create_list_params(
+            params, param_lists, param_cols, **kwargs
+        )
+        return cls(list_params, **kwargs)
 
     @staticmethod
     def create_list_params(
         params: Dict[str, MFParam],
         param_lists: list,
+        param_cols: list,
         **kwargs,
     ) -> Dict[str, MFParam]:
         """Create the param dictionary"""
         idx = 0
         list_params = dict()
         for param_name, param in params.items():
-            if type(param) is MFDouble:
+            if param_name == "aux" or param_name == "boundname":
+                continue
+            shape = None
+            if param_cols[idx] == 1:
+                shape = len(param_lists[idx])
+            else:
+                shape = (len(param_lists[idx]), param_cols[idx])
+            if type(param) is MFArray and param.type == "double":
                 list_params[param_name] = MFArray(
-                    shape=len(param_lists[idx]),
+                    shape=shape,
                     array=np.array(param_lists[idx], dtype=np.float64),
                     how=MFArrayType.internal,
                     factor=1.0,
                     path=None,
                     **kwargs,
                 )
-            elif type(param) is MFInteger:
+            elif type(param) is MFArray and param.type == "integer":
                 list_params[param_name] = MFArray(
-                    shape=len(param_lists[idx]),
+                    shape=shape,
                     array=np.array(param_lists[idx], dtype=np.int32),
                     how=MFArrayType.internal,
                     factor=1,
@@ -406,7 +453,7 @@ def create_list_params(
             else:
                 list_params[param_name] = MFScalarList(
                     value=param_lists[idx],
-                    type=type(param),
+                    # type=type(param),
                     **kwargs,
                 )
 
@@ -427,5 +474,9 @@ def write(self, f, **kwargs):
         for i in range(count):
             line = f"{PAD}"
             for name, param in self.params.items():
-                line += f"{param.value[i]}\t"
+                if isinstance(param.value[i], np.ndarray):
+                    for v in param.value[i]:
+                        line += f"{v}\t"
+                else:
+                    line += f"{param.value[i]}\t"
             f.write(line + "\n")

flopy4/ispec/exg_gwfgwf.py

@@ -463,7 +463,7 @@ class ExgGwfgwf(MFPackage):
     )
 
     aux = MFArray(
-        type = "array",
+        type = "double",
         block = "exchangedata",
         shape = "(naux)",
         reader = "urword",