fix(GwfGwfExchangeModule): correct specific discharge calculation for…

… LGR (#56)

* refactor(pre-commit.py): use OrderdedDict to load and write code.json

Updated the code.json writer to use an ordered dictionary so that the order is always the same.  This will minimize the number of changes that appear in the changelog due to ordering changes that do not change content.

* fix(GwfGwfExchangeModule): specific discharge not updated correctly for LGR

The sign for the face normal needed to be flipped for model 2 so that specific discharge would be calculated correctly.  Added a test to make sure that the velocities are correct.  Addresses #51.

* fix(GwfGwfExchangeModule): specific discharge not updated correctly for LGR

Updated the release notes and added some docstrings to the autotest.

autotest/test_gwf_npf04_spdis.py

@@ -0,0 +1,235 @@
+"""
+MODFLOW 6 Autotest
+Test the specific discharge calculation for an LGR-like simulation that has
+a parent model and a child model.  The child model is inset into the parent
+model, but they both have the same resolution, so it is essentially a simple
+3D grid.  The child qx velocity should be the same as the qx velocity in
+the parent grid.  The heads are also compared.
+
+"""
+
+import os
+import sys
+import numpy as np
+
+try:
+    import flopy
+    from flopy.utils.lgrutil import Lgr
+except:
+    msg = 'Error. FloPy package is not available.\n'
+    msg += 'Try installing using the following command:\n'
+    msg += ' pip install flopy'
+    raise Exception(msg)
+
+from framework import testing_framework
+from simulation import Simulation
+
+ex = ['npf04']
+exdirs = []
+for s in ex:
+    exdirs.append(os.path.join('temp', s))
+
+namea = 'a'
+nameb = 'b'
+
+
+def get_model(idx, dir):
+
+    # grid properties
+    nlay = 3
+    nrow = 6
+    ncol = 6
+    delr = 100.
+    delc = 100.
+    top = 300.
+    botm = [200., 100., 0.]
+
+    # hydraulic properties
+    hk = 1.
+    vk = 1.
+
+    # Set the idomain of the parent model in order to
+    # define where the child model will be located
+    idomain = np.ones((nlay, nrow, ncol), dtype=int)
+    idomain[:, 2:4, 2:4] = 0
+
+    ncpp = 1
+    ncppl = [1, 1, 1]
+    lgr = Lgr(nlay, nrow, ncol, delr, delc, top, botm,
+              idomain, ncpp, ncppl)
+
+    name = ex[idx]
+
+    # build MODFLOW 6 files
+    # create simulation
+    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+                                 exe_name='mf6',
+                                 sim_ws=dir)
+
+    # create tdis package
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    # create gwf model
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=namea, save_flows=True)
+
+    # create iterative model solution and register the gwf model with it
+    ims = flopy.mf6.ModflowIms(sim, outer_hclose=1e-9, inner_hclose=1.e-9)
+
+    # dis
+    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
+                                  delr=delr, delc=delc,
+                                  top=top, botm=botm, idomain=idomain)
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, pname='ic', strt=top)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True,
+                                  icelltype=0, k=hk, k33=vk)
+
+    # chd
+    chdspd = []
+    for k in range(nlay):
+        for i in range(nrow):
+            chdspd.append([(k, i, 0), 1.])
+            chdspd.append([(k, i, ncol - 1), 6.])
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(gwf, pname='oc',
+                                budget_filerecord='{}.cbc'.format(namea),
+                                head_filerecord='{}.hds'.format(namea),
+                                headprintrecord=[('COLUMNS', 10, 'WIDTH', 15,
+                                                  'DIGITS', 6, 'GENERAL')],
+                                saverecord=[('HEAD', 'ALL'),
+                                            ('BUDGET', 'ALL')],
+                                printrecord=[('HEAD', 'ALL'),
+                                             ('BUDGET', 'ALL')])
+
+    # create child gwf model
+    cgwf = flopy.mf6.ModflowGwf(sim, modelname=nameb, save_flows=True)
+    cnlay, cnrow, cncol = lgr.get_shape()
+    cdelr, cdelc = lgr.get_delr_delc()
+    ctop, cbotm = lgr.get_top_botm()
+    xorigin, yorigin = lgr.get_lower_left()
+    cidomain = lgr.get_idomain()
+    cdis = flopy.mf6.ModflowGwfdis(cgwf, nlay=cnlay, nrow=cnrow, ncol=cncol,
+                                   delr=cdelr, delc=cdelc,
+                                   top=ctop, botm=cbotm, idomain=cidomain,
+                                   xorigin=xorigin, yorigin=yorigin)
+    cic = flopy.mf6.ModflowGwfic(cgwf, pname='ic', strt=top)
+    cnpf = flopy.mf6.ModflowGwfnpf(cgwf, save_specific_discharge=True,
+                                   icelltype=0, k=hk, k33=vk)
+    oc = flopy.mf6.ModflowGwfoc(cgwf, pname='oc',
+                                budget_filerecord='{}.cbc'.format(nameb),
+                                head_filerecord='{}.hds'.format(nameb),
+                                headprintrecord=[('COLUMNS', 10, 'WIDTH', 15,
+                                                  'DIGITS', 6, 'GENERAL')],
+                                saverecord=[('HEAD', 'ALL'),
+                                            ('BUDGET', 'ALL')],
+                                printrecord=[('HEAD', 'ALL'),
+                                             ('BUDGET', 'ALL')])
+
+    # exchange information
+    exchangedata = lgr.get_exchange_data(angldegx=True, cdist=True)
+    nexg = len(exchangedata)
+    gwfe = flopy.mf6.ModflowGwfgwf(sim, exgtype='gwf6-gwf6',
+                                   exgmnamea='a', exgmnameb='b',
+                                   nexg=nexg,
+                                   auxiliary=[('angldegx', 'cdist')],
+                                   exchangedata=exchangedata)
+
+    return sim
+
+
+def build_models():
+    for idx, dir in enumerate(exdirs):
+        sim = get_model(idx, dir)
+        sim.write_simulation()
+    return
+
+
+def qxqyqz(fname, nlay, nrow, ncol):
+    nodes = nlay * nrow * ncol
+    cbb = flopy.utils.CellBudgetFile(fname, precision='double')
+    spdis = cbb.get_data(text='DATA-SPDIS')[0]
+    qx = np.ones((nodes), dtype=np.float) * 1.e30
+    qy = np.ones((nodes), dtype=np.float) * 1.e30
+    qz = np.ones((nodes), dtype=np.float) * 1.e30
+    n0 = spdis['node'] - 1
+    qx[n0] = spdis['qx']
+    qy[n0] = spdis['qy']
+    qz[n0] = spdis['qz']
+    qx = qx.reshape(nlay, nrow, ncol)
+    qy = qy.reshape(nlay, nrow, ncol)
+    qz = qz.reshape(nlay, nrow, ncol)
+    qx = np.ma.masked_equal(qx, 1.e30)
+    qy = np.ma.masked_equal(qy, 1.e30)
+    qz = np.ma.masked_equal(qz, 1.e30)
+    return qx, qy, qz
+
+
+def eval_mf6(sim):
+    print('evaluating head and qx in parent and child models...')
+
+    # make sure parent head is same as child head in same column
+    fname = os.path.join(sim.simpath, '{}.hds'.format(namea))
+    hdobj = flopy.utils.HeadFile(fname)
+    ha = hdobj.get_data()
+    fname = os.path.join(sim.simpath, '{}.hds'.format(nameb))
+    hdobj = flopy.utils.HeadFile(fname)
+    hb = hdobj.get_data()
+    msg = 'Heads should be the same {} {}'.format(ha[0, 1, 2], hb[0, 0, 0])
+    assert np.allclose(ha[0, 1, 2], hb[0, 0, 0]), msg
+
+    # make sure specific discharge is calculated correctly for child and
+    # parent models (even though child model has same resolution as parent
+    fname = os.path.join(sim.simpath, '{}.cbc'.format(namea))
+    nlaya, nrowa, ncola = ha.shape
+    qxa, qya, qza = qxqyqz(fname, nlaya, nrowa, ncola)
+    fname = os.path.join(sim.simpath, '{}.cbc'.format(nameb))
+    nlayb, nrowb, ncolb = hb.shape
+    qxb, qyb, qzb = qxqyqz(fname, nlayb, nrowb, ncolb)
+    msg = 'qx should be the same {} {}'.format(qxa[0, 2, 1], qxb[0, 0, 0])
+    assert np.allclose(qxa[0, 2, 1], qxb[0, 0, 0]), msg
+
+    return
+
+
+# - No need to change any code below
+def test_mf6model():
+    # initialize testing framework
+    test = testing_framework()
+
+    # build the models
+    build_models()
+
+    # run the test models
+    for idx, dir in enumerate(exdirs):
+        yield test.run_mf6, Simulation(dir, exfunc=eval_mf6, idxsim=idx)
+
+    return
+
+
+def main():
+    # initialize testing framework
+    test = testing_framework()
+
+    # build the models
+    build_models()
+
+    # run the test models
+    for idx, dir in enumerate(exdirs):
+        sim = Simulation(dir, exfunc=eval_mf6, idxsim=idx)
+        test.run_mf6(sim)
+
+    return
+
+
+
+if __name__ == "__main__":
+    # print message
+    print('standalone run of {}'.format(os.path.basename(__file__)))
+
+    # run main routine
+    main()

doc/ReleaseNotes/ReleaseNotes.tex

@@ -123,7 +123,8 @@ \section{History}
 \item Addressed issue with pointing contiguous pointer vectors/arrays to non-contiguous pointer vectors/arrays that caused code compilation failure with gfortran-8. A consequence of addressing this issue is that all pointer vectors/arrays that are allocated or pointed to using the memory manager must be defined to be contiguous.
 \item Corrected a problem with the reading of grid data from a binary file, in which the program was reading a binary header for each row of data.
 \item Added a new error check for very small time steps.  If the value of starttime plus the time step length is equal to endtime, then the time step is too small to be differentiated by the program based on the precision of floating point numbers.  The program will terminate with an error in this case.  The program will also terminate if the storage package and a transient stress period has a time step length of zero.
-\item The observation package was modified to use non-advancing output instead of fixed length strings when writing ascii output. The previous use of fixed length strings resulted in truncation of ascii observation output when the product of user-specified \texttt{digits} + 7 and the number of observations exceeded 5000. 
+\item The observation package was modified to use non-advancing output instead of fixed length strings when writing ascii output. The previous use of fixed length strings resulted in truncation of ascii observation output when the product of user-specified \texttt{digits} + 7 and the number of observations exceeded 5000.
+\item Corrected an error in the GWF-GWF Exchange module that caused the specific discharge values in the child model to be calculated incorrectly.  The calculation was incorrect because the face normal for the child model was pointing toward the center of the cell instead of outward.
 \end{itemize}
 
 \underline{STRESS PACKAGES}