Fixes: #603, #606, #614, and #616

autotest/t504_test.py

@@ -318,6 +318,16 @@ def test006_gwf3():
     # load simulation
     sim = MFSimulation.load(model_name, 'mf6', exe_name, pth)
 
+    model = sim.get_model()
+    disu = model.get_package('disu')
+    # test switching disu array to internal array
+    disu.ja = disu.ja.array
+    # test writing hwva and cl12 arrays out to different locations
+    disu.hwva = {'filename': 'flow.disu.hwva_new.dat', 'factor': 1.0,
+               'data': disu.hwva.array}
+    disu.cl12 = {'filename': 'flow.disu.cl12_new.dat', 'factor': 1.0,
+               'data': disu.cl12.array}
+
     # make temp folder to save simulation
     sim.simulation_data.mfpath.set_sim_path(run_folder)
     # write simulation to new location
@@ -858,12 +868,12 @@ def test027_timeseriestest():
 
 
 if __name__ == '__main__':
+    test006_gwf3()
     test001a_tharmonic()
     test001e_uzf_3lay()
     test003_gwfs_disv()
     test005_advgw_tidal()
     test006_2models_mvr()
-    test006_gwf3()
     test027_timeseriestest()
     test036_twrihfb()
     test045_lake1ss_table()

flopy/export/utils.py

@@ -600,7 +600,7 @@ def transient2d_export(f, t2d, **kwargs):
 
     if isinstance(f, str) and f.lower().endswith(".shp"):
         array_dict = {}
-        for kper in range(t2d.model.modeltime.sim_time.nper):
+        for kper in range(t2d.model.modeltime.nper):
             u2d = t2d[kper]
             name = '{}_{}'.format(
                 shapefile_utils.shape_attr_name(u2d.name), kper + 1)

flopy/mf6/coordinates/modelgrid.py

@@ -433,7 +433,7 @@ def grid_type_consistent(self):
 
     def get_connections_array(self):
         if self.grid_type() == DiscretizationType.DISU:
-            return np.arange(1, self.num_connections() + 1, 1, np.int)
+            return np.arange(1, self.num_connections() + 1, 1, np.int32)
         else:
             except_str = 'ERROR: Can not get connections arrays for model ' \
                          '"{}" Only DISU (unstructured) grids ' \
@@ -443,10 +443,10 @@ def get_connections_array(self):
 
     def get_horizontal_cross_section_dim_arrays(self):
         if self.grid_type() == DiscretizationType.DIS:
-            return [np.arange(1, self.num_rows() + 1, 1, np.int),
-                    np.arange(1, self.num_columns() + 1, 1, np.int)]
+            return [np.arange(1, self.num_rows() + 1, 1, np.int32),
+                    np.arange(1, self.num_columns() + 1, 1, np.int32)]
         elif self.grid_type() == DiscretizationType.DISV:
-            return [np.arange(1, self.num_cells_per_layer() + 1, 1, np.int)]
+            return [np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)]
         elif self.grid_type() == DiscretizationType.DISU:
             except_str = 'ERROR: Can not get horizontal plane arrays for ' \
                          'model "{}" DISU grid.  DISU grids do not support ' \
@@ -464,23 +464,23 @@ def get_model_dim(self):
 
     def get_model_dim_arrays(self):
         if self.grid_type() == DiscretizationType.DIS:
-            return [np.arange(1, self.num_layers() + 1, 1, np.int),
-                    np.arange(1, self.num_rows() + 1, 1, np.int),
-                    np.arange(1, self.num_columns() + 1, 1, np.int)]
+            return [np.arange(1, self.num_layers() + 1, 1, np.int32),
+                    np.arange(1, self.num_rows() + 1, 1, np.int32),
+                    np.arange(1, self.num_columns() + 1, 1, np.int32)]
         elif self.grid_type() == DiscretizationType.DISV:
-            return [np.arange(1, self.num_layers() + 1, 1, np.int),
-                    np.arange(1, self.num_cells_per_layer() + 1, 1, np.int)]
+            return [np.arange(1, self.num_layers() + 1, 1, np.int32),
+                    np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)]
         elif self.grid_type() == DiscretizationType.DISU:
-            return [np.arange(1, self.num_cells() + 1, 1, np.int)]
+            return [np.arange(1, self.num_cells() + 1, 1, np.int32)]
 
     def get_row_array(self):
-        return np.arange(1, self.num_rows() + 1, 1, np.int)
+        return np.arange(1, self.num_rows() + 1, 1, np.int32)
 
     def get_column_array(self):
-        return np.arange(1, self.num_columns() + 1, 1, np.int)
+        return np.arange(1, self.num_columns() + 1, 1, np.int32)
 
     def get_layer_array(self):
-        return np.arange(1, self.num_layers() + 1, 1, np.int)
+        return np.arange(1, self.num_layers() + 1, 1, np.int32)
 
     def get_horizontal_cross_section_dim_names(self):
         if self.grid_type() == DiscretizationType.DIS:

flopy/mf6/data/dfn/gwf-disu.dfn

@@ -111,6 +111,7 @@ reader readarray
 longname grid connectivity
 description is a list of cell number (n) followed by its connecting cell numbers (m) for each of the m cells connected to cell n. The number of values to provide for cell n is IAC(n).  This list is sequentially provided for the first to the last cell. The first value in the list must be cell n itself, and the remaining cells must be listed in an increasing order (sorted from lowest number to highest).  Note that the cell and its connections are only supplied for the GWF cells and their connections to the other GWF cells.  Also note that the JA list input may be divided such that every node and its connectivity list can be on a separate line for ease in readability of the file. To further ease readability of the file, the node number of the cell whose connectivity is subsequently listed, may be expressed as a negative number, the sign of which is subsequently converted to positive by the code.
 numeric_index true
+jagged_array iac
 
 block connectiondata
 name ihc
@@ -119,6 +120,7 @@ shape (nja)
 reader readarray
 longname connection type
 description is an index array indicating the direction between node n and all of its m connections.  If IHC = 0 then cell n and cell m are connected in the vertical direction.  Cell n overlies cell m if the cell number for n is less than m; cell m overlies cell n if the cell number for m is less than n.  If IHC = 1 then cell n and cell m are connected in the horizontal direction.  If IHC = 2 then cell n and cell m are connected in the horizontal direction, and the connection is vertically staggered.  A vertically staggered connection is one in which a cell is horizontally connected to more than one cell in a horizontal connection.
+jagged_array iac
 
 block connectiondata
 name cl12
@@ -127,6 +129,7 @@ shape (nja)
 reader readarray
 longname connection lengths
 description is the array containing connection lengths between the center of cell n and the shared face with each adjacent m cell.
+jagged_array iac
 
 block connectiondata
 name hwva
@@ -135,6 +138,7 @@ shape (nja)
 reader readarray
 longname connection lengths
 description is a symmetric array of size NJA.  For horizontal connections, entries in HWVA are the horizontal width perpendicular to flow.  For vertical connections, entries in HWVA are the vertical area for flow.  Thus, values in the HWVA array contain dimensions of both length and area.  Entries in the HWVA array have a one-to-one correspondence with the connections specified in the JA array.  Likewise, there is a one-to-one correspondence between entries in the HWVA array and entries in the IHC array, which specifies the connection type (horizontal or vertical).  Entries in the HWVA array must be symmetric; the program will terminate with an error if the value for HWVA for an n to m connection does not equal the value for HWVA for the corresponding n to m connection.
+jagged_array iac
 
 block connectiondata
 name angldegx
@@ -144,6 +148,7 @@ shape (nja)
 reader readarray
 longname angle of face normal to connection
 description is the angle (in degrees) between the horizontal x-axis and the outward normal to the face between a cell and its connecting cells. The angle varies between zero and 360.0 degrees, where zero degrees points in the positive x-axis direction, and 90 degrees points in the positive y-axis direction.  ANGLDEGX is only needed if horizontal anisotropy is specified in the NPF Package, if the XT3D option is used in the NPF Package, or if the SAVE\_SPECIFIC\_DISCHARGE option is specifed in the NPF Package.  ANGLDEGX does not need to be specified if these conditions are not met.  ANGLDEGX is of size NJA; values specified for vertical connections and for the diagonal position are not used.  Note that ANGLDEGX is read in degrees, which is different from MODFLOW-USG, which reads a similar variable (ANGLEX) in radians.
+jagged_array iac
 
 # --------------------- gwf disu vertices ---------------------
 

flopy/mf6/data/mfdata.py

@@ -448,7 +448,10 @@ def _get_internal_formatting_string(self, layer):
         if storage.data_structure_type != DataStructureType.recarray:
             int_format.append('FACTOR')
             if layer_storage.factor is not None:
-                int_format.append(str(layer_storage.factor))
+                if data_type == DatumType.integer:
+                    int_format.append(str(int(layer_storage.factor)))
+                else:
+                    int_format.append(str(layer_storage.factor))
             else:
                 if data_type == DatumType.double_precision:
                     int_format.append('1.0')
@@ -475,8 +478,12 @@ def _get_external_formatting_string(self, layer, ext_file_action):
         ext_format = ['OPEN/CLOSE', "'{}'".format(ext_file_path)]
         if storage.data_structure_type != DataStructureType.recarray:
             if layer_storage.factor is not None:
+                data_type = self.structure.get_datum_type(return_enum_type=True)
                 ext_format.append('FACTOR')
-                ext_format.append(str(layer_storage.factor))
+                if data_type == DatumType.integer:
+                    ext_format.append(str(int(layer_storage.factor)))
+                else:
+                    ext_format.append(str(layer_storage.factor))
         if layer_storage.binary:
             ext_format.append('(BINARY)')
         if layer_storage.iprn is not None:

flopy/mf6/data/mfdataarray.py

@@ -565,6 +565,10 @@ def get_file_entry(self, layer=None,
                                       inspect.stack()[0][3], type_,
                                       value_, traceback_, None,
                                       self._simulation_data.debug, ex)
+            if self.structure.data_item_structures[0].numeric_index or \
+                    self.structure.data_item_structures[0].is_cellid:
+                # for cellid and numeric indices convert from 0 base to 1 based
+                data = abs(data) + 1
             file_entry_array.append('{}{}{}{}\n'.format(indent,
                                                         self.structure.name,
                                                         indent,
@@ -680,9 +684,8 @@ def _get_file_entry_layer(self, layer, data_indent, storage_type,
             # internal data header + data
             format_str = self._get_internal_formatting_string(layer).upper()
             lay_str = self._get_data_layer_string(layer, data_indent).upper()
-            file_entry = '{}{}{}\n{}{}'.format(file_entry, indent_string,
-                                               format_str, indent_string,
-                                               lay_str)
+            file_entry = '{}{}{}\n{}'.format(file_entry, indent_string,
+                                               format_str, lay_str)
         elif storage_type == DataStorageType.internal_constant:
             #  constant data
             try:
@@ -728,8 +731,6 @@ def _get_file_entry_layer(self, layer, data_indent, storage_type,
         return file_entry
 
     def _get_data_layer_string(self, layer, data_indent):
-        layer_data_string = ['']
-        line_data_count = 0
         # iterate through data layer
         try:
             data = self._get_storage_obj().get_data(layer, False)
@@ -744,42 +745,10 @@ def _get_data_layer_string(self, layer, data_indent):
                                   inspect.stack()[0][3], type_,
                                   value_, traceback_, comment,
                                   self._simulation_data.debug, ex)
-        data_iter = datautil.PyListUtil.next_item(data)
-        indent_str = self._simulation_data.indent_string
-        for item, last_item, new_list, nesting_change in data_iter:
-            # increment data/layer counts
-            line_data_count += 1
-            try:
-                data_lyr = to_string(item, self._data_type,
-                                     self._simulation_data,
-                                     self._data_dimensions)
-            except Exception as ex:
-                type_, value_, traceback_ = sys.exc_info()
-                comment = 'Could not convert data "{}" of type "{}" to a ' \
-                          'string.'.format(item, self._data_type)
-                raise MFDataException(self.structure.get_model(),
-                                      self.structure.get_package(),
-                                      self._path,
-                                      'converting data',
-                                      self.structure.name,
-                                      inspect.stack()[0][3], type_,
-                                      value_, traceback_, comment,
-                                      self._simulation_data.debug, ex)
-            layer_data_string[-1] = '{}{}{}'.format(layer_data_string[-1],
-                                                    indent_str,
-                                                    data_lyr)
-            if self._simulation_data.wrap_multidim_arrays and \
-                    (line_data_count == self._simulation_data.
-                        max_columns_of_data or last_item):
-                layer_data_string.append('{}'.format(data_indent))
-                line_data_count = 0
-        if len(layer_data_string) > 0:
-            # clean up the text at the end of the array
-            layer_data_string[-1] = layer_data_string[-1].strip()
-        if len(layer_data_string) == 1:
-            return '{}{}\n'.format(data_indent, layer_data_string[0].rstrip())
-        else:
-            return '\n'.join(layer_data_string)
+        file_access = MFFileAccessArray(self.structure, self._data_dimensions,
+                                        self._simulation_data, self._path,
+                                        self._current_key)
+        return file_access.get_data_string(data, self._data_type, data_indent)
 
     def _resolve_layer_index(self, layer, allow_multiple_layers=False):
         # handle layered vs non-layered data

flopy/mf6/data/mfdatalist.py

@@ -178,7 +178,7 @@ def to_array(self, kper=0, mask=False):
                 raise Exception("MfList: something bad happened")
 
         for name, arr in arrays.items():
-            cnt = np.zeros(shape, dtype=np.float)
+            cnt = np.zeros(shape, dtype=np.float64)
             #print(name,kper)
             for sp_rec in sarr:
                 if sp_rec is not None:

flopy/mf6/data/mfdatascalar.py

@@ -80,17 +80,17 @@ def plotable(self):
     @property
     def dtype(self):
         if self.structure.type == DatumType.double_precision:
-            return np.float32
+            return np.float64
         elif self.structure.type == DatumType.integer:
-            return np.int
+            return np.int32
         elif self.structure.type == DatumType.recarray or \
                 self.structure.type == DatumType.record or \
                 self.structure.type == DatumType.repeating_record:
             for data_item_struct in self.structure.data_item_structures:
                 if data_item_struct.type == DatumType.double_precision:
-                    return np.float32
+                    return np.float64
                 elif data_item_struct.type == DatumType.integer:
-                    return np.int
+                    return np.int32
         return None
 
     def has_data(self):
@@ -169,7 +169,7 @@ def set_data(self, data):
 
     def add_one(self):
         datum_type = self.structure.get_datum_type()
-        if datum_type == int or datum_type == np.int:
+        if datum_type == int or datum_type == np.int32:
             if self._get_storage_obj().get_data() is None:
                 try:
                     self._get_storage_obj().set_data(1)

flopy/mf6/data/mfdatastorage.py

@@ -1092,11 +1092,12 @@ def store_external(self, file_path, layer=None, multiplier=None,
                         self.data_dimensions.structure, self.data_dimensions,
                         self._simulation_data, self._data_path,
                         self._stress_period)
+                    str_layered = self.data_dimensions.structure.layered
                     file_access.write_binary_file(
                         data, fp, text, self._model_or_sim.modeldiscrit,
                         self._model_or_sim.modeltime,
                         stress_period=self._stress_period, precision='double',
-                        write_multi_layer=(layer is None))
+                        write_multi_layer=(layer is None and str_layered))
                 else:
                     file_access = MFFileAccessArray(
                         self.data_dimensions.structure, self.data_dimensions,
@@ -1573,7 +1574,7 @@ def _build_full_data(self, apply_multiplier=False):
         if dimensions[0] < 0:
             return None
         all_none = True
-        np_data_type =  self.data_dimensions.structure.get_datum_type()
+        np_data_type = self.data_dimensions.structure.get_datum_type()
         full_data = np.full(dimensions, np.nan,
                             self.data_dimensions.structure.get_datum_type(True))
         is_aux = self.data_dimensions.structure.name == 'aux'
@@ -1675,7 +1676,9 @@ def _fill_const_layer(self, layer):
         if data_dimensions[0] < 0:
             return ls.data_const_value
         else:
-            return np.full(data_dimensions, ls.data_const_value[0])
+            data_type = self.data_dimensions.structure. \
+                get_datum_type(numpy_type=True)
+            return np.full(data_dimensions, ls.data_const_value[0], data_type)
 
     def _is_type(self, data_item, data_type):
         if data_type == DatumType.string or data_type == DatumType.keyword:

flopy/mf6/data/mfdatautil.py

@@ -3,7 +3,7 @@
 from copy import deepcopy
 from ..mfbase import MFDataException, FlopyException
 from .mfstructure import DatumType
-from ...utils.datautil import PyListUtil
+from ...utils.datautil import PyListUtil, DatumUtil
 import struct
 
 
@@ -122,6 +122,8 @@ def to_string(val, data_type, sim_data, data_dim, is_cellid=False,
             else:
                 return sim_data.sci_format_str.format(val)
     elif is_cellid or (possible_cellid and isinstance(val, tuple)):
+        if DatumUtil.is_int(val):
+            return str(val + 1)
         if len(val) > 0 and val[0] == 'none':
             # handle case that cellid is 'none'
             return val[0]
@@ -526,7 +528,7 @@ def _build_layer(self, data_type, data_storage_type, default_value,
                                    data_type)
         elif data_storage_type == DataStorageType.internal_constant:
             if default_value is None:
-                if data_type == np.int:
+                if data_type == np.int32:
                     data = 0
                 else:
                     data = 0.0