feat(maw): program head limit and rate scaling to work for injection …

…wells

Updated code, added a new test, updated readme and mf6io files.

Also fix a few minor typos in the definition files.  And remove the PRECISION keyword from observation options since it is not supported.

README.md

@@ -7,7 +7,7 @@
 
 ## Automated Testing Status on Travis-CI
 
-### Version 6.0.3 develop — build 38
+### Version 6.0.3 develop — build 56
 [![Build Status](https://travis-ci.org/MODFLOW-USGS/modflow6.svg?branch=develop)](https://travis-ci.org/MODFLOW-USGS/modflow6)
 
 ## Introduction
@@ -31,7 +31,7 @@ MODFLOW 6 is the latest core version of MODFLOW. It synthesizes many of the capa
 
 #### ***Software/Code citation for MODFLOW 6:***
 
-[Langevin, C.D., Hughes, J.D., Banta, E.R., Provost, A.M., Niswonger, R.G., and Panday, Sorab, 2019, MODFLOW 6 Modular Hydrologic Model version 6.0.3 — develop: U.S. Geological Survey Software Release, 07 January 2019, https://doi.org/10.5066/F76Q1VQV](https://doi.org/10.5066/F76Q1VQV)
+[Langevin, C.D., Hughes, J.D., Banta, E.R., Provost, A.M., Niswonger, R.G., and Panday, Sorab, 2019, MODFLOW 6 Modular Hydrologic Model version 6.0.3 — develop: U.S. Geological Survey Software Release, 07 February 2019, https://doi.org/10.5066/F76Q1VQV](https://doi.org/10.5066/F76Q1VQV)
 
 ## Instructions for building definition files for new packages
 

autotest/test_gwf_maw01.py

@@ -36,10 +36,8 @@ def get_model(idx, dir):
     perlen = [1., 1., 1.]
     nstp = [1, 1, 1]
     tsmult = [1., 1., 1.]
-    steady = [True, True, True]
     lenx = 300.
     delr = delc = lenx / float(nrow)
-    botm = [0.]
     strt = 100.
     hnoflo = 1e30
     hdry = -1e30

autotest/test_gwf_maw02.py

@@ -36,10 +36,10 @@
 shape3d = (nlay, nrow, ncol)
 size3d = nlay * nrow * ncol
 nper = 5
-perlen = [1. for n in range(nper)]
-nstp = [1 for n in range(nper)]
-tsmult = [1. for n in range(nper)]
-steady = [True for n in range(nper)]
+perlen = nper * [1.]
+nstp = nper * [1]
+tsmult = nper * [1.]
+steady = nper * [True]
 lenx = 300.
 delr = delc = lenx / float(nrow)
 botm = [0.]

autotest/test_gwf_maw03.py

@@ -0,0 +1,246 @@
+"""
+MODFLOW 6 Autotest
+Test the MAW HEAD_LIMIT and RATE_SCALING options for injection wells.  These
+options were not originally supported in MODFLOW 6.  They were added for
+version 6.0.4.
+
+"""
+
+import os
+import sys
+import numpy as np
+
+try:
+    import pymake
+except:
+    msg = 'Error. Pymake package is not available.\n'
+    msg += 'Try installing using the following command:\n'
+    msg += ' pip install https://github.com/modflowpy/pymake/zipball/master'
+    raise Exception(msg)
+
+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)
+
+from framework import testing_framework
+from simulation import Simulation
+
+ex = ['maw03a', 'maw03b', 'maw03c']
+exdirs = []
+for s in ex:
+    exdirs.append(os.path.join('temp', s))
+
+# maw settings for runs a, b, and c
+mawsetting_a = [(0, 'rate', 2000.), ]
+mawsetting_b = [(0, 'rate', 2000.), (0, 'head_limit', 0.4)]
+mawsetting_c = [(0, 'rate', 2000.), (0, 'rate_scaling', 0.0, 1.0)]
+mawsettings = [mawsetting_a, mawsetting_b, mawsetting_c]
+
+def get_model(idx, dir):
+
+    nlay, nrow, ncol = 1, 101, 101
+    nper = 1
+    perlen = [1000.]
+    nstp = [50]
+    tsmult = [1.2]
+    delr = delc = 142.
+    top = 0.
+    botm = [-1000.]
+    strt = 0.
+    hk = 10.
+
+    nouter, ninner = 100, 100
+    hclose, rclose, relax = 1e-6, 1e-6, 1.
+
+    tdis_rc = []
+    for i in range(nper):
+        tdis_rc.append((perlen[i], nstp[i], tsmult[i]))
+
+    name = ex[idx]
+
+    # build MODFLOW 6 files
+    ws = dir
+    sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws)
+
+    # create tdis package
+    tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
+                                 nper=nper, perioddata=tdis_rc)
+
+    # create gwf model
+    gwf = flopy.mf6.MFModel(sim, model_type='gwf6', modelname=name,
+                            model_nam_file='{}.nam'.format(name))
+
+    # create iterative model solution and register the gwf model with it
+    ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
+                               outer_hclose=hclose,
+                               outer_maximum=nouter,
+                               under_relaxation='NONE',
+                               inner_maximum=ninner,
+                               inner_hclose=hclose, rcloserecord=rclose,
+                               linear_acceleration='CG',
+                               scaling_method='NONE',
+                               reordering_method='NONE',
+                               relaxation_factor=relax)
+    sim.register_ims_package(ims, [gwf.name])
+
+    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
+                                  delr=delr, delc=delc,
+                                  top=top, botm=botm,
+                                  idomain=1,
+                                  fname='{}.dis'.format(name))
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt,
+                                fname='{}.ic'.format(name))
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True,
+                                  icelltype=1,
+                                  k=hk,
+                                  k33=hk,
+                                  fname='{}.npf'.format(name))
+
+    # storage
+    sto = flopy.mf6.ModflowGwfsto(gwf, save_flows=True,
+                                  iconvert=0,
+                                  ss=1.e-5, sy=0.1,
+                                  steady_state={0: False},
+                                  transient={0: True},
+                                  fname='{}.sto'.format(name))
+
+    # MAW
+    opth = '{}.maw.obs'.format(name)
+    wellbottom = -1000
+    wellrecarray = [[0, 0.15, wellbottom, 0., 'THIEM', 1]]
+    wellconnectionsrecarray = [[0, 0, (0, 50, 50), 0., wellbottom, 0., 0.0]]
+    wellperiodrecarray = mawsettings[idx]
+    maw = flopy.mf6.ModflowGwfmaw(gwf, fname='{}.maw'.format(name),
+                                  print_input=True, print_head=True,
+                                  print_flows=True, save_flows=True,
+                                  obs_filerecord=opth,
+                                  packagedata=wellrecarray,
+                                  connectiondata=wellconnectionsrecarray,
+                                  perioddata=wellperiodrecarray)
+    mawo_dict = {}
+    mawo_dict['{}.maw.obs.csv'.format(name)] = [('m1head', 'head', (0,)),
+                                                ('m1rate', 'rate', (0,))]       #is this index one-based? Not if in a tuple
+    maw_obs = flopy.mf6.ModflowUtlobs(gwf,
+                                      fname=opth,
+                                      parent_file=maw, digits=15,
+                                      print_input=True,
+                                      continuous=mawo_dict)
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(gwf,
+                                budget_filerecord='{}.cbc'.format(name),
+                                head_filerecord='{}.hds'.format(name),
+                                headprintrecord=[
+                                    ('COLUMNS', ncol, 'WIDTH', 15,
+                                     'DIGITS', 6, 'GENERAL')],
+                                saverecord=[('HEAD', 'ALL')],
+                                printrecord=[('HEAD', 'ALL'),
+                                             ('BUDGET', 'ALL')],
+                                fname='{}.oc'.format(name))
+
+    # head observations
+    obs_data0 = [('head_well_cell', 'HEAD', (0, 0, 0))]
+    obs_recarray = {'{}.obs.csv'.format(name): obs_data0}
+    obs = flopy.mf6.ModflowUtlobs(gwf, pname='head_obs',
+                                  fname='{}.obs'.format(name),
+                                  digits=15, print_input=True,
+                                  continuous=obs_recarray)
+
+
+    return sim
+
+
+def build_models():
+    for idx, dir in enumerate(exdirs):
+        sim = get_model(idx, dir)
+        sim.write_simulation()
+    return
+
+
+def eval_maw(sim):
+    print('evaluating MAW heads...')
+
+    # MODFLOW 6 maw results
+    idx = sim.idxsim
+    name = ex[idx]
+    fpth = os.path.join(sim.simpath, '{}.maw.obs.csv'.format(name))
+    try:
+        tc = np.genfromtxt(fpth, names=True, delimiter=',')
+    except:
+        assert False, 'could not load data from "{}"'.format(fpth)
+
+
+    if idx == 0:
+
+        # M1RATE should be 2000.
+        msg = 'The injection rate should be 2000. for all times'
+        assert tc['M1RATE'].min() == tc['M1RATE'].max() == 2000, msg
+
+    elif idx == 1:
+
+        # M1RATE should have a minimum value less than 200 and
+        # M1HEAD should not exceed 0.400001
+        msg = ('Injection rate should fall below 200 and the head should not'
+              'exceed 0.4')
+        assert tc['M1RATE'].min() < 200., msg
+        assert tc['M1HEAD'].max() < 0.400001, msg
+
+    elif idx == 2:
+
+        # M1RATE should have a minimum value less than 800
+        # M1HEAD should not exceed 1.0.
+        msg = ('Min injection rate should be less than 800 and well '
+               'head should not exceed 1.0')
+        assert tc['M1RATE'].min() < 800. and tc['M1HEAD'].max() < 1., msg
+
+    else:
+
+        assert False, 'Test error.  idx {} not being tested.'.format(idx)
+
+    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_maw, 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_maw, 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()

code.json

@@ -1,38 +1,38 @@
 [
     {
-        "status": "Release Candidate",
+        "status": "Release Candidate", 
         "languages": [
             "Fortran2008"
-        ],
-        "repositoryURL": "https://code.usgs.gov/usgs/modflow/modflow6.git",
-        "disclaimerURL": "https://code.usgs.gov/usgs/modflow/modflow6/blob/master/DISCLAIMER.md",
+        ], 
+        "repositoryURL": "https://code.usgs.gov/usgs/modflow/modflow6.git", 
+        "disclaimerURL": "https://code.usgs.gov/usgs/modflow/modflow6/blob/master/DISCLAIMER.md", 
         "tags": [
-            "MODFLOW",
+            "MODFLOW", 
             "groundwater model"
-        ],
-        "vcs": "git",
-        "name": "modflow6",
-        "downloadURL": "https://code.usgs.gov/usgs/modflow/modflow6/archive/master.zip",
+        ], 
+        "vcs": "git", 
+        "name": "modflow6", 
+        "downloadURL": "https://code.usgs.gov/usgs/modflow/modflow6/archive/master.zip", 
         "contact": {
-            "name": "Christian D. Langevin",
+            "name": "Christian D. Langevin", 
             "email": "langevin@usgs.gov"
-        },
-        "laborHours": -1,
-        "version": "6.0.3.38",
+        }, 
+        "laborHours": -1, 
+        "version": "6.0.3.56", 
         "date": {
-            "metadataLastUpdated": "2019-01-07"
-        },
-        "organization": "U.S. Geological Survey",
+            "metadataLastUpdated": "2019-02-07"
+        }, 
+        "organization": "U.S. Geological Survey", 
         "permissions": {
             "licenses": [
                 {
-                    "URL": "https://code.usgs.gov/usgs/modflow/modflow6/blob/master/LICENSE.md",
+                    "URL": "https://code.usgs.gov/usgs/modflow/modflow6/blob/master/LICENSE.md", 
                     "name": "Public Domain, CC0-1.0"
                 }
-            ],
+            ], 
             "usageType": "openSource"
-        },
-        "homepageURL": "https://code.usgs.gov/usgs/modflow/modflow6/",
+        }, 
+        "homepageURL": "https://code.usgs.gov/usgs/modflow/modflow6/", 
         "description": "MODFLOW is the USGS's modular hydrologic model. MODFLOW is considered an international standard for simulating and predicting groundwater conditions and groundwater/surface-water interactions."
     }
 ]

doc/ReleaseNotes/ReleaseNotes.tex

@@ -122,7 +122,7 @@ \section{History}
 \begin{itemize}
 \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 Added a new error check for very small time steps.  If the value of the starting time is equal to the ending time (starting time plus the time step length), 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 with 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 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}
@@ -136,7 +136,7 @@ \section{History}
 
 \underline{ADVANCED STRESS PACKAGES}
 \begin{itemize}
-\item -
+\item Modified the Multi-Aquifer Well (MAW) Package so that the HEAD\_LIMIT and RATE\_SCALING options work for injection wells.  Prior to this change, these options only worked for extraction wells.  These options can be used to reduce or even shut off well injection as the head in the well rises above user-specified levels.
 \item -
 \item - 
 \end{itemize}

doc/mf6io/gwf/binaryoutput.tex

@@ -580,7 +580,7 @@ \subsection{Observation Output File}
 
 \begin{description} \itemsep0pt \parskip0pt \parsep0pt
 \item \texttt{TYPE} (bytes 1--4 of Record 1) is ``cont `` ---  ``cont'' indicates the file contains continuous observations ;
-\item \texttt{PRECISION} (bytes 6--11 of Record 1) is either ``single'' or ``double'' --- ``single'' indicates that floating-point values are written in single precision (4 bytes), and ``double'' indicates double precision (8 bytes);
+\item \texttt{PRECISION} (bytes 6--11 of Record 1) will always be ``double'' to indicate that floating-point values are written in double precision (8 bytes);
 \item \texttt{LENOBSNAME} (bytes 12--15 of Record 1) is an integer indicating the number of characters used to store each observation name in following records (in the initial release of MODFLOW~6, LENOBSNAME equals 40);
 \item \texttt{NOBS} (4-byte integer) is the number of observations recorded in the file;
 \item \texttt{OBSNAME} (LENOBSNAME bytes) is an observation name;

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

@@ -378,7 +378,7 @@ tagged false
 in_record true
 reader urword
 longname lake connection type
-description character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include:  VERTICAL--character keyword to indicate the lake-GWF connection is vertical  and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKCONN = 1) and there must be a lake table associated with each embedded lake.
+description character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include:  VERTICAL--character keyword to indicate the lake-GWF connection is vertical  and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKECONN = 1) and there must be a lake table associated with each embedded lake.
 
 block connectiondata
 name bedleak