feat(disu idomain): add support for disu idomain (#261)

* feat(disu): adding idomain capability to disu

* couple of fixes to get disu working for non-reduced grids

* Found another place in disu where a conversion to user node number is required

* change variable name

todo list:
still need to consider converting iac to ia right away
change usr to input for variable names so as not to conflict with iausr/jausr in connections
add more comprehensive test
add idomain validity check
support idomain < 0 values

* more cleanup on disu idomain support

* more disu cleanup and idomain testing

* cleanup comments

* Updated release notes

autotest/test_gwf_disu01.py

@@ -0,0 +1,140 @@
+"""
+MODFLOW 6 Autotest
+Test to make sure that disu is working correctly
+
+"""
+
+import os
+import shutil
+import subprocess
+import numpy as np
+from nose.tools import raises
+
+try:
+    import flopy
+except:
+    msg = 'Error. FloPy package is not available.\n'
+    msg += 'Try installing using the following command:\n'
+    msg += ' pip install flopy'
+    raise Exception(msg)
+
+import targets
+
+mf6_exe = os.path.abspath(targets.target_dict['mf6'])
+testname = 'gwf_disu01'
+testdir = os.path.join('temp', testname)
+if not os.path.isdir(testdir):
+    os.mkdir(testdir)
+everything_was_successful = True
+
+
+def run_mf6(argv, ws):
+    buff = []
+    proc = subprocess.Popen(argv,
+                            stdout=subprocess.PIPE,
+                            stderr=subprocess.PIPE,
+                            cwd=ws)
+    result, error = proc.communicate()
+    if result is not None:
+        c = result.decode('utf-8')
+        c = c.rstrip('\r\n')
+        print('{}'.format(c))
+        buff.append(c)
+
+    return proc.returncode, buff
+
+
+def test_disu_simple():
+    from disu_util import get_disu_kwargs
+    name = 'disu01a'
+    ws = os.path.join(testdir, name)
+    nlay = 3
+    nrow = 3
+    ncol = 3
+    delr = 10. * np.ones(ncol)
+    delc = 10. * np.ones(nrow)
+    top = 0
+    botm = [-10, -20, -30]
+    disukwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm)
+    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+                                 exe_name=mf6_exe, sim_ws=ws)
+    tdis = flopy.mf6.ModflowTdis(sim)
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
+    ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY')
+    disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs)
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=0.)
+    npf = flopy.mf6.ModflowGwfnpf(gwf)
+    spd = {0:[[(0, ), 1.], [(nrow * ncol - 1, ), 0.]]}
+    chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd)
+    sim.write_simulation()
+    sim.run_simulation()
+    return
+
+
+def test_disu_idomain_simple():
+    from disu_util import get_disu_kwargs
+    name = 'disu01b'
+    ws = os.path.join(testdir, name)
+    nlay = 3
+    nrow = 3
+    ncol = 3
+    delr = 10. * np.ones(ncol)
+    delc = 10. * np.ones(nrow)
+    top = 0
+    botm = [-10, -20, -30]
+    idomain = np.ones(nlay*nrow*ncol, dtype=np.int)
+    idomain[1] = 0
+    disukwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm)
+    disukwargs['idomain'] = idomain
+    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+                                 exe_name=mf6_exe, sim_ws=ws)
+    tdis = flopy.mf6.ModflowTdis(sim)
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
+    ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY')
+    disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs)
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=0.)
+    npf = flopy.mf6.ModflowGwfnpf(gwf)
+    spd = {0:[[(0, ), 1.], [(nrow * ncol - 1, ), 0.]]}
+    chd = flopy.mf6.modflow.ModflowGwfchd(gwf, stress_period_data=spd)
+    oc = flopy.mf6.modflow.ModflowGwfoc(gwf,
+                                        budget_filerecord='{}.bud'.format(name),
+                                        head_filerecord='{}.hds'.format(name),
+                                        saverecord=[('HEAD', 'LAST'),
+                                                    ('BUDGET', 'LAST')],)
+    sim.write_simulation()
+    sim.run_simulation()
+
+    # check binary grid file
+    fname = os.path.join(ws, name + '.disu.grb')
+    grbobj = flopy.utils.MfGrdFile(fname)
+    nodes = grbobj._datadict['NODES']
+    ia = grbobj._datadict['IA']
+    ja = grbobj._datadict['JA']
+    assert nodes == disukwargs['nodes']
+    assert np.array_equal(ia[0: 4], np.array([1, 4, 4, 7]))
+    assert np.array_equal(ja[:6], np.array([ 1, 4, 10, 3, 6, 12]))
+    assert ia[-1] == 127
+    assert ia.shape[0] == 28, 'ia should have size of 28'
+    assert ja.shape[0] == 126, 'ja should have size of 126'
+
+    # load head array and ensure nodata value in second cell
+    fname = os.path.join(ws, name + '.hds')
+    hdsobj = flopy.utils.HeadFile(fname)
+    head = hdsobj.get_alldata().flatten()
+    assert head[1] == 1.e30
+
+    # load flowja to make sure it is the right size
+    fname = os.path.join(ws, name + '.bud')
+    budobj = flopy.utils.CellBudgetFile(fname, precision='double')
+    flowja = budobj.get_data(text='FLOW-JA-FACE')[0].flatten()
+    assert flowja.shape[0] == 126
+
+    return
+
+
+if __name__ == "__main__":
+    # print message
+    print('standalone run of {}'.format(os.path.basename(__file__)))
+
+    test_disu_simple()
+    test_disu_idomain_simple()

doc/ReleaseNotes/ReleaseNotes.tex

@@ -120,22 +120,23 @@ \section{History}
 \underline{BASIC FUNCTIONALITY}
 \begin{itemize}
 \item Added an error check to the DISU Package that ensures that an underlying cell has a top elevation that is less than or equal to the bottom of an overlying cell.  An underlying cell is one in which the IHC value for the connection is zero and the connecting node number is greater than the cell node number.
+\item Added restricted IDOMAIN support for DISU grids.  Users can specify an optional IDOMAIN in the DISU Package input file.  IDOMAIN values must be zero or one.  Vertical pass-through cells (specified with an IDOMAIN value of -1 in the DIS or DISV Package input files) are not supported for DISU.   
 \item NPF Package will now write a message to the GWF Model list file to indicate when the SAVE\_SPECIFIC\_DISCHARGE option is invoked.
 \item Added two new options to the NPF Package.  The K22OVERK option allows the user to enter the anisotropy ratio for K22.  If activated, the K22 values entered by the user in the NPF input file will be multiplied by the K values entered in the NPF input file.  The K33OVERK option allows the user to enter the anisotropy ratio for K33.  If activated, the K33 values entered by the user in the NPF input file will be multiplied by the K values entered in the NPF input file.  With this K33OVERK option, for example, the user can specify a value of 0.1 for K33 and have all K33 values be one tenth of the values specified for K.  The program will terminate with an error if these options are invoked, but arrays for K22 and/or K33 are not provided in the NPF input file.
 \end{itemize}
 
 \underline{STRESS PACKAGES}
 \begin{itemize}
 \item There was an error in the calculation of the segmented evapotranspiration rate for the case where the rate did not decrease with depth.  There was another error in which PETM0 was being used as the evapotranspiration rate at the surface instead of the proportion of the evapotranspiration rate at the surface.
-\item Corrected several error messages issued by the SFR Package that were not formatted correctly.  
-\item Fixed a bug in which the lake stage stable would sometimes result in touching numbers.  This only occurred for negative lake stages.
 \end{itemize}
 
 \underline{ADVANCED STRESS PACKAGES}
 \begin{itemize}
 \item Corrected the way auxiliary variables are handled for the advanced packages.  In some cases, values for auxiliary variables were not being correctly written to the GWF Model budget file or to the advanced package budget file.  A consistent approach for updating and saving auxiliary variables was implemented for the MAW, SFR, LAK, and UZF Packages.
 \item The user guide was updated to include a missing laksetting that was omitted from the PERIOD block.  The laksetting description now includes an INFLOW option; a description for INFLOW is also now included.
-\item For the advanced stress packages, values assigned to the auxiliary variables were not written correctly to the GWF Model budget file, but the values were correct in the advanced package budget file.  Program was modified so that auxiliary variables written to the GWF Model budget file are correct.
+\item For the advanced stress packages, values assigned to the auxiliary variables were not written correctly to the GWF Model budget file, but the values were correct in the advanced package budget file.  Program was modified so that auxiliary variables are correctly written to the GWF Model budget file.
+\item Corrected several error messages issued by the SFR Package that were not formatted correctly.  
+\item Fixed a bug in which the lake stage stable would sometimes result in touching numbers.  This only occurred for negative lake stages.
 \end{itemize}
 
 \underline{SOLUTION}

doc/mf6io/mf6ivar/dfn/gwf-disu.dfn

@@ -92,6 +92,15 @@ reader readarray
 longname cell surface area
 description is the cell surface area (in plan view).
 
+block griddata
+name idomain
+type integer
+shape (nodes)
+reader readarray
+layered true
+optional true
+longname idomain existence array
+description is an optional array that characterizes the existence status of a cell.  If the IDOMAIN array is not specified, then all model cells exist within the solution.  If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation.  Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution.  If the IDOMAIN value for a cell is 1, the cell exists in the simulation.  IDOMAIN values of -1 cannot be specified for the DISU Package.
 
 # --------------------- gwf disu connectiondata ---------------------