write runnable baseline gwf simulations (#43)

* write runnable gwf base simulation

* remove temp path checks

* add gwf disu test

* add gwf disv test

* misc cleanup

---------

Co-authored-by: mjreno <mjreno@IGSAAA071L01144.gs.doi.net>

Author: mjreno

docs/examples/array_example.py

@@ -62,11 +62,12 @@
 constant = data_path / "constant.txt"
 external = data_path / "external.txt"
 shape = (1000, 100)
+dtype = "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=dtype, header=False)
 
 # Get values
 
@@ -87,7 +88,7 @@
 plt.colorbar()
 
 fhandle = open(constant)
-cmfa = MFArray.load(fhandle, data_path, shape, header=False)
+cmfa = MFArray.load(fhandle, data_path, shape, type=dtype, header=False)
 cvals = cmfa.value
 plt.imshow(cvals[0:100])
 plt.colorbar()
@@ -110,7 +111,7 @@
 # External
 
 fhandle = open(external)
-emfa = MFArray.load(fhandle, data_path, shape, header=False)
+emfa = MFArray.load(fhandle, data_path, shape, type=dtype, header=False)
 evals = emfa.value
 evals
 
@@ -135,7 +136,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=dtype, header=False, layered=True
+)
 vals = ilmfa.value
 
 ilmfa._value  # internal storage
@@ -182,7 +185,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=dtype, header=False, layered=True
+)
 
 clmfa._value
 
@@ -235,7 +240,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=dtype, 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])
@@ -446,20 +456,30 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
                     nlay = params.get("dimensions").get("nlay")
                     nrow = params.get("dimensions").get("nrow")
                     ncol = params.get("dimensions").get("ncol")
+                    shape = (nlay, nrow, ncol)
+                elif "disv" in mempath.split("/"):
+                    nlay = params.get("dimensions").get("nlay")
                     ncpl = params.get("dimensions").get("ncpl")
-                    nodes = params.get("dimensions").get("nodes")
-                    if nrow and ncol:
-                        shape = (nlay, nrow, ncol)
-                    elif ncpl:
+                    nvert = params.get("dimensions").get("nvert")
+                    if shape == "(ncpl)":
+                        shape = ncpl
+                    elif shape == "(ncpl, nlay)":
                         shape = (nlay, ncpl)
-                    elif nodes:
+                    elif shape == "(nvert)":
+                        shape = nvert
+                elif "disu" in mempath.split("/"):
+                    nodes = params.get("dimensions").get("nodes")
+                    nja = params.get("dimensions").get("nja")
+                    if "nodes" in shape:
                         shape = nodes
+                    elif "nja" in shape:
+                        shape = nja
         if layered:
             nlay = shape[0]
             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 +494,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 +517,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 +553,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.
@@ -532,11 +564,11 @@ def read_array(f):
             pos = f.tell()
             line = f.readline()
             line = line_strip(line)
-            if not re.match("^[0-9. ]+$", line):
+            if not re.match("^[-0-9. ]+$", line):
                 f.seek(pos, 0)
                 break
             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
 
@@ -235,15 +236,20 @@ def load(cls, f, **kwargs):
             line = f.readline()
             if line == "":
                 raise ValueError("Early EOF, aborting")
-            if line == "\n":
+            if line == "\n" or line.lstrip().startswith("#"):
                 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)