Transport (#63)

* feat(GwtSrbModule): created new sorption package for transport model

First crack at the sorption package.  The GwtSrbModule has been implemented, but has not yet been added to the gwt model.  The source compiles and existing tests run and pass so the new code doesn't break anything.

* feat(GwtSrbModule) : finished up much of the Sorption Package structure

Sorption package is up and running and has a test, input instructions, and flopy support.  Right now it only has equilibrium-control linear sorption.  Should also add Freundlich, Langmuir and the nonequilibrium options that are in MT3D.  And what about dual-domain?  Does that go here?  And there may still be a question about the volume that should be used for calculating sorption.  Right now, we are using Vcell, but this might need another look.  Also still need to add and test observations.

* test(GwtSrbModule): removed a hard coded assert False

Removed a statement that was left in by accident.

* docs(readme): fix directory name

* fix(GwfModule): fix to apply ptc to problems with the STO package

Fix to apply ptc during steady-state stress periods to models using
the Newton-Raphson formulation and including the STO package.
Previously ptc was only applied to models using the Newton-Raphson
formulation and no STO package - this is not the intended behaviour.
Also migrated DEV_NO_PTC option to NO_PTC option. Added autotest
(test_gwf_ptc01) that tests that the same results are achieved for
problems with and without the STO packages (this test is identical to
MODFLOW-NWT Problem 3 - low recharge).

closes #47

* 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.

* refactor(InputOutputModule): refactor dclosetest (#49)

Refactor dclosetest and replace with IS_SAME in GenericUtilities.
Convert IS_SAME to logical function and allow passing of evaluation
value (eps). If eps is not passed then DSAME is used. Modified calls
to dclosetest to IS_SAME but did not pass an eps value (unlike what
was done previously). All tests pass so current tests are not sensitive
to a passed eps value. Need to monitor this for timeseries functionality
which used dclosetest.

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

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.

* perf(GwfModule): improvements to ptc (#53)

Small improvements to ptc to reduce the ptcdel value loaded on the
diagonal.

* fix(ObsModule): implemented non-advancing output (#55)

Modified the pbservation process 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 digits + 7 and the number of
observations exceeded 5000.

Closes #54

* 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.

* docs(dfn): fix 3 dfn typos (#59)

- `double` → `double precision`
- `valid_values` → `valid`
- `in_record = false` → `in_record false`

* docs(mf6ivar): readme.md updated to include valid keyword (#60)

* 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.

* docs(mf6ivar): readme.md needed a description of the valid keyword

* fix(ssm): added trap for case where flows are not available to transport model

CONTRIBUTING.md

@@ -131,7 +131,7 @@ from the main (upstream) repository:
 ## <a name="rules"></a> Coding Rules
 To ensure consistency throughout the source code, keep these rules in mind as you are working:
 
-* All features or bug fixes **must be tested** by one or more specs (unit-tests).
+* All features or bug fixes **must be tested** by one or more specs (unit-tests and/or integration/regression-tests).
 
 ## <a name="commit"></a> Commit Message Guidelines
 
@@ -175,7 +175,7 @@ If the commit reverts a previous commit, it should begin with `revert: `, follow
 ### Type
 Must be one of the following:
 
-* **ci**: Changes to our CI configuration files and scripts (example scopes: Travis, Circle, BrowserStack, SauceLabs)
+* **ci**: Changes to our CI configuration files and scripts (example scopes: Travis)
 * **docs**: Documentation only changes
 * **feat**: A new feature
 * **fix**: A bug fix
@@ -185,13 +185,10 @@ Must be one of the following:
 * **test**: Adding missing tests or correcting existing tests
 
 ### Scope
-The scope should be the name of the MODFLOW 6 module/class affected (as perceived by the person reading the changelog generated from commit messages.
+The scope should be the name of the MODFLOW 6 module/class affected (as perceived by the person reading the changelog generated from commit messages).
 
 There are currently a few exceptions to the "use module/class name" rule:
 
-* **packaging**: used for changes that change the npm package layout in all of our packages, e.g.
-  public path changes, package.json changes done to all packages, d.ts file/format changes, changes
-  to bundles, etc.
 * **releasenotes**: used for updating the release notes
 * **readme**: used for updating the release notes in README.md
 * **changelog**: used for updating the release notes in CHANGELOG.md

DEVELOPER.md

@@ -81,7 +81,7 @@ To run tests:
 
 ```shell
 # Go to the autotest directory
-cd modflow6/autotests
+cd modflow6/autotest
 
 # Run all modflow6 tests (including building executables and the mfio documentation - requires installation of LaTeX
 nosetests -v                 

autotest/simulation.py

@@ -66,6 +66,7 @@ def __init__(self, name, exfunc=None, exe_dict=None, htol=None,
         self.simpath = None
         self.inpt = None
         self.outp = None
+        self.coutp = None
 
         # set htol for comparisons
         if htol is None:
@@ -95,13 +96,16 @@ def set_model(self, pth):
 
         self.simpath = pth
 
-        # get MNODFLOW6 output file names
+        # get MODFLOW 6 output file names
         fpth = os.path.join(pth, 'mfsim.nam')
         mf6inp, mf6outp = pymake.get_mf6_files(fpth)
         self.outp = mf6outp
 
         # determine comparison model
         self.action = pymake.get_mf6_comparison(pth)
+        if self.action is not None:
+            if 'mf6' in self.action:
+                cinp, self.coutp = pymake.get_mf6_files(fpth)
 
     def setup(self, src, dst):
         msg = sfmt.format('Setup test', self.name)
@@ -169,8 +173,11 @@ def run(self):
                     cpth = os.path.join(self.simpath, self.action)
                     key = self.action.lower().replace('.cmp', '')
                     exe = os.path.abspath(targets.target_dict[key])
-                    npth = pymake.get_namefiles(cpth)[0]
-                    nam = os.path.basename(npth)
+                    if 'mf6' in key:
+                        nam = None
+                    else:
+                        npth = pymake.get_namefiles(cpth)[0]
+                        nam = os.path.basename(npth)
                     self.nam_cmp = nam
                     try:
                         success_cmp, buff = flopy.run_model(exe, nam,
@@ -249,7 +256,17 @@ def compare(self):
                                 os.path.basename(pth))
                             print(txt)
                     else:
-                        files2.append(None)
+                        if self.coutp is not None:
+                            for file2 in self.coutp:
+                                ext = os.path.splitext(file2)[1][1:]
+
+                                if ext.lower() in ['hds', 'hed', 'bhd', 'ahd']:
+                                    # simulation file
+                                    pth = os.path.join(cpth, file2)
+                                    files2.append(pth)
+
+                        else:
+                            files2.append(None)
 
             if self.nam_cmp is None:
                 pth = None

autotest/test_gwf_maw01.py

@@ -22,7 +22,7 @@
 from simulation import Simulation
 
 ex = ['maw01', 'maw01nwt', 'maw01nwtur']
-newtonoptions = [None, [''], ['UNDER_RELAXATION']]
+newtonoptions = [None, [''], ['NEWTON', 'UNDER_RELAXATION']]
 exdirs = []
 for s in ex:
     exdirs.append(os.path.join('temp', s))

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()