fix(data check): Minimum number of data columns check fixed (#1062)

autotest/t505_test.py

@@ -37,7 +37,8 @@
 
 try:
     import shapefile
-    if int(shapefile.__version__.split('.')[0]) < 2:
+
+    if int(shapefile.__version__.split(".")[0]) < 2:
         shapefile = None
 except ImportError:
     shapefile = None
@@ -1134,6 +1135,28 @@ def test005_advgw_tidal():
         timeseries=ts_dict,
     )
 
+    # test nseg = 1
+    evt_period = ModflowGwfevt.stress_period_data.empty(model, 150, nseg=1)
+    for col in range(0, 10):
+        for row in range(0, 15):
+            evt_period[0][col * 15 + row] = (
+                (0, row, col),
+                50.0,
+                0.0004,
+                10.0,
+                None,
+            )
+    evt_package_test = ModflowGwfevt(
+        model,
+        print_input=True,
+        print_flows=True,
+        save_flows=True,
+        maxbound=150,
+        nseg=1,
+        stress_period_data=evt_period,
+    )
+    evt_package_test.remove()
+
     # test empty
     evt_period = ModflowGwfevt.stress_period_data.empty(model, 150, nseg=3)
     for col in range(0, 10):
@@ -1264,7 +1287,15 @@ def test005_advgw_tidal():
             ("rv2-upper", "RIV", "riv2_upper"),
             ("rv-2-7-4", "RIV", (0, 6, 3)),
             ("rv2-8-5", "RIV", (0, 6, 4)),
-            ("rv-2-9-6", "RIV", (0, 5, 5,)),
+            (
+                "rv-2-9-6",
+                "RIV",
+                (
+                    0,
+                    5,
+                    5,
+                ),
+            ),
         ],
         "riv_flowsA.csv": [
             ("riv1-3-1", "RIV", (0, 2, 0)),
@@ -2816,6 +2847,196 @@ def test028_sfr():
     return
 
 
+def test_transport():
+    # init paths
+    test_ex_name = "test_transport"
+    name = "mst03"
+
+    pth = os.path.join(
+        "..", "examples", "data", "mf6", "create_tests", test_ex_name
+    )
+    run_folder = os.path.join(cpth, test_ex_name)
+    if not os.path.isdir(run_folder):
+        os.makedirs(run_folder)
+
+    expected_output_folder = os.path.join(pth, "expected_output")
+    expected_head_file = os.path.join(expected_output_folder, "gwf_mst03.hds")
+    expected_conc_file = os.path.join(expected_output_folder, "gwt_mst03.unc")
+
+    ws = os.path.join('temp', 't505', test_ex_name)
+    laytyp = [1]
+    ss = [1.e-10]
+    sy = [0.1]
+    nlay, nrow, ncol = 1, 1, 1
+
+    nper = 2
+    perlen = [2., 2.]
+    nstp = [14, 14]
+    tsmult = [1., 1.]
+    delr = 10.
+    delc = 10.
+    top = 10.
+    botm = [0.]
+    strt = top
+    hk = 1.0
+
+    nouter, ninner = 100, 300
+    hclose, rclose, relax = 1e-6, 1e-6, 0.97
+
+    tdis_rc = []
+    for idx in range(nper):
+        tdis_rc.append((perlen[idx], nstp[idx], tsmult[idx]))
+    idx = 0
+
+    # build MODFLOW 6 files
+    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+                                 exe_name='mf6',
+                                 sim_ws=ws)
+    # create tdis package
+    tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
+                                 nper=nper, perioddata=tdis_rc)
+
+    # create gwf model
+    gwfname = 'gwf_' + name
+    newtonoptions = ['NEWTON', 'UNDER_RELAXATION']
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname,
+                               newtonoptions=newtonoptions,)
+
+    # create iterative model solution and register the gwf model with it
+    imsgwf = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
+                                  outer_dvclose=hclose,
+                                  outer_maximum=nouter,
+                                  under_relaxation='DBD',
+                                  under_relaxation_theta=0.7,
+                                  inner_maximum=ninner,
+                                  inner_dvclose=hclose, rcloserecord=rclose,
+                                  linear_acceleration='BICGSTAB',
+                                  scaling_method='NONE',
+                                  reordering_method='NONE',
+                                  relaxation_factor=relax,
+                                  filename='{}.ims'.format(gwfname))
+    sim.register_ims_package(imsgwf, [gwf.name])
+
+    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
+                                  delr=delr, delc=delc,
+                                  top=top, botm=botm,
+                                  idomain=np.ones((nlay, nrow, ncol),
+                                                  dtype=int))
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=False,
+                                  icelltype=laytyp[idx],
+                                  k=hk, k33=hk)
+    # storage
+    sto = flopy.mf6.ModflowGwfsto(gwf, save_flows=False,
+                                  iconvert=laytyp[idx],
+                                  ss=ss[idx], sy=sy[idx],
+                                  steady_state={0: False},
+                                  transient={0: True})
+
+    # wel files
+    welspdict = {0: [[(0, 0, 0), -25., 0.]],
+                 1: [[(0, 0, 0), 25., 0.]]}
+    wel = flopy.mf6.ModflowGwfwel(gwf, print_input=True, print_flows=True,
+                                  stress_period_data=welspdict,
+                                  save_flows=False,
+                                  auxiliary='CONCENTRATION', pname='WEL-1')
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(gwf,
+                                budget_filerecord='{}.cbc'.format(gwfname),
+                                head_filerecord='{}.hds'.format(gwfname),
+                                headprintrecord=[
+                                    ('COLUMNS', 10, 'WIDTH', 15,
+                                     'DIGITS', 6, 'GENERAL')],
+                                saverecord=[('HEAD', 'ALL')],
+                                printrecord=[('HEAD', 'ALL'),
+                                             ('BUDGET', 'ALL')])
+
+    # create gwt model
+    gwtname = 'gwt_' + name
+    gwt = flopy.mf6.MFModel(sim, model_type='gwt6', modelname=gwtname,
+                            model_nam_file='{}.nam'.format(gwtname))
+
+    # create iterative model solution and register the gwt model with it
+    imsgwt = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
+                                  outer_dvclose=hclose,
+                                  outer_maximum=nouter,
+                                  under_relaxation='NONE',
+                                  inner_maximum=ninner,
+                                  inner_dvclose=hclose, rcloserecord=rclose,
+                                  linear_acceleration='BICGSTAB',
+                                  scaling_method='NONE',
+                                  reordering_method='NONE',
+                                  relaxation_factor=relax,
+                                  filename='{}.ims'.format(gwtname))
+    sim.register_ims_package(imsgwt, [gwt.name])
+
+    dis = flopy.mf6.ModflowGwtdis(gwt, nlay=nlay, nrow=nrow, ncol=ncol,
+                                  delr=delr, delc=delc,
+                                  top=top, botm=botm,
+                                  idomain=1,
+                                  filename='{}.dis'.format(gwtname))
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwtic(gwt, strt=100.)
+
+    # advection
+    adv = flopy.mf6.ModflowGwtadv(gwt, scheme='UPSTREAM',
+                                filename='{}.adv'.format(gwtname))
+
+    # mass storage and transfer
+    mst = flopy.mf6.ModflowGwtmst(gwt, porosity=sy[idx],
+                                filename='{}.mst'.format(gwtname))
+
+    # sources
+    sourcerecarray = [('WEL-1', 'AUX', 'CONCENTRATION')]
+    ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray,
+                                filename='{}.ssm'.format(gwtname))
+
+    # output control
+    oc = flopy.mf6.ModflowGwtoc(gwt,
+                                budget_filerecord='{}.cbc'.format(gwtname),
+                                concentration_filerecord='{}.ucn'.format(gwtname),
+                                concentrationprintrecord=[
+                                    ('COLUMNS', 10, 'WIDTH', 15,
+                                     'DIGITS', 6, 'GENERAL')],
+                                saverecord=[('CONCENTRATION', 'ALL')],
+                                printrecord=[('CONCENTRATION', 'ALL'),
+                                             ('BUDGET', 'ALL')])
+
+    # GWF GWT exchange
+    gwfgwt = flopy.mf6.ModflowGwfgwt(sim, exgtype='GWF6-GWT6',
+                                     exgmnamea=gwfname, exgmnameb=gwtname,
+                                     filename='{}.gwfgwt'.format(name))
+
+    # write MODFLOW 6 files
+    sim.write_simulation()
+
+    # run simulation
+    if run:
+        sim.run_simulation()
+
+        # compare output to expected results
+        head_file = os.path.join(os.getcwd(), expected_head_file)
+        head_new = os.path.join(run_folder, "gwf_mst03.hds")
+        outfile = os.path.join(run_folder, "head_compare.dat")
+        assert pymake.compare_heads(
+            None, None, files1=head_file, files2=head_new, outfile=outfile
+        )
+        conc_file = os.path.join(os.getcwd(), expected_conc_file)
+        conc_new = os.path.join(run_folder, "gwt_mst03.ucn")
+        assert pymake.compare_concs(
+            None, None, files1=conc_file, files2=conc_new, outfile=outfile
+        )
+
+        # clean up
+        sim.delete_output_files()
+
+
 if __name__ == "__main__":
     np001()
     np002()
@@ -2827,3 +3048,4 @@ def test028_sfr():
     test028_sfr()
     test035_fhb()
     test050_circle_island()
+    test_transport()

flopy/mf6/coordinates/modeldimensions.py

@@ -4,6 +4,7 @@
 
 """
 
+import sys
 from .simulationtime import SimulationTime
 from .modelgrid import UnstructuredModelGrid, ModelGrid
 from ..mfbase import StructException, FlopyException, VerbosityLevel
@@ -81,6 +82,7 @@ def get_data_shape(
         data=None,
         data_item_num=None,
         repeating_key=None,
+        min_size=False,
     ):
         return self.get_model_dim(data_item_num).get_data_shape(
             self.structure,
@@ -89,6 +91,7 @@ def get_data_shape(
             data,
             self.package_dim.package_path,
             repeating_key=repeating_key,
+            min_size=min_size,
         )
 
     def model_subspace_size(self, subspace_string="", data_item_num=None):
@@ -403,6 +406,7 @@ def get_data_shape(
         path=None,
         deconstruct_axis=True,
         repeating_key=None,
+        min_size=False,
     ):
         if structure is None:
             raise FlopyException(
@@ -492,6 +496,7 @@ def get_data_shape(
                 path,
                 deconstruct_axis,
                 repeating_key=repeating_key,
+                min_size=min_size,
             )
             if self.locked and consistent_shape:
                 self.stored_shapes[data_item.path] = (
@@ -509,6 +514,7 @@ def _resolve_data_item_shape(
         path=None,
         deconstruct_axis=True,
         repeating_key=None,
+        min_size=False,
     ):
         if isinstance(data, tuple):
             data = [data]
@@ -560,7 +566,7 @@ def _resolve_data_item_shape(
                         result = self.resolve_exp(
                             item,
                             self._find_in_dataset(
-                                data_set_struct, item[0], data
+                                data_set_struct, item[0], data, min_size
                             ),
                         )
                         if result:
@@ -584,7 +590,8 @@ def _resolve_data_item_shape(
                         elif DatumUtil.is_int(item[0]):
                             shape_dimensions.append(int(item[0]))
                         else:
-                            # try to resolve dimension within the existing block
+                            # try to resolve dimension within the
+                            # existing block
                             result = self.simulation_data.mfdata.find_in_path(
                                 parent_path, item[0]
                             )
@@ -684,7 +691,7 @@ def resolve_exp(self, expression, value):
             return value
 
     @staticmethod
-    def _find_in_dataset(data_set_struct, item, data):
+    def _find_in_dataset(data_set_struct, item, data, min_size=False):
         if data is not None:
             # find the current data item in data_set_struct
             for index, data_item in zip(
@@ -695,12 +702,22 @@ def _find_in_dataset(data_set_struct, item, data):
                     data_item.name.lower() == item.lower()
                     and len(data[0]) > index
                 ):
-                    # always use the maximum value
-                    max_val = 0
-                    for data_line in data:
-                        if data_line[index] > max_val:
-                            max_val = data_line[index]
-                    return max_val
+                    if min_size:
+                        # use the minimum value
+                        min_val = sys.maxsize
+                        for data_line in data:
+                            if data_line[index] < min_val:
+                                min_val = data_line[index]
+                        if min_val == sys.maxsize:
+                            return 0
+                        return min_val
+                    else:
+                        # use the maximum value
+                        max_val = 0
+                        for data_line in data:
+                            if data_line[index] > max_val:
+                                max_val = data_line[index]
+                        return max_val
         return None
 
     @staticmethod