fix(ObsModule): implemented non-advancing output

README.md

@@ -7,8 +7,8 @@
 
 ## Automated Testing Status on Travis-CI
 
-### Version 6.0.3 feat-ptcdist — build 30
-[![Build Status](https://travis-ci.org/MODFLOW-USGS/modflow6.svg?branch=feat-ptcdist)](https://travis-ci.org/MODFLOW-USGS/modflow6)
+### Version 6.0.3 obs-outlinelen — build 31
+[![Build Status](https://travis-ci.org/MODFLOW-USGS/modflow6.svg?branch=obs-outlinelen)](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, 2018, MODFLOW 6 Modular Hydrologic Model version 6.0.3 — feat-ptcdist: U.S. Geological Survey Software Release, 29 October 2018, 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, 2018, MODFLOW 6 Modular Hydrologic Model version 6.0.3 — obs-outlinelen: U.S. Geological Survey Software Release, 30 October 2018, https://doi.org/10.5066/F76Q1VQV](https://doi.org/10.5066/F76Q1VQV)
 
 
 ## Instructions for building definition files for new packages

autotest/test_gwf_utl03_obs01.py

@@ -0,0 +1,244 @@
+import os
+import sys
+import numpy as np
+
+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 = ['utl03_obs']
+exdirs = []
+for s in ex:
+    exdirs.append(os.path.join('temp', s))
+ddir = 'data'
+
+# temporal discretization
+nper = 2
+tdis_rc = [(1., 1, 1.),
+           (1., 1, 1.)]
+
+# spatial discretization data
+nlay, nrow, ncol = 1, 21, 21
+top = 200.
+botm = 0.
+
+# hydraulic properties
+laytyp = 1
+hk = 1.
+delr = delc = 100.
+strt = botm + 20.
+
+# constant head data
+chdlocl = [(0, i, 0) for i in range(nrow)]
+chdlocr = [(0, i, ncol - 1) for i in range(nrow)]
+chdl = 100.
+chdr = 75.
+
+c60 = []
+c61 = []
+for loc in chdlocl:
+    c60.append([loc, chdl])
+    c61.append([loc, chdl + 10])
+for loc in chdlocr:
+    c60.append([loc, chdr])
+    c61.append([loc, chdr - 10])
+cd6 = {0: c60, 1: c61}
+
+# gwf obs
+obs_data0 = [('h{:04d}'.format(i + 1), 'HEAD', (0, 10, 10)) for i in
+             range(1000)]
+obs_data1 = [('h{:04d}'.format(i + 1001), 'HEAD', (0, 1, 1)) for i in
+             range(737)]
+
+# solver data
+nouter, ninner = 100, 300
+hclose, rclose, relax = 1e-6, 0.01, 1.
+
+
+def build_mf6(idx, ws, binaryobs=True):
+    name = ex[idx]
+
+    # build MODFLOW 6 files
+    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+                                 exe_name='mf6',
+                                 sim_ws=ws)
+    # create tdis package
+    flopy.mf6.ModflowTdis(sim, time_units='DAYS',
+                          nper=nper, perioddata=tdis_rc)
+
+    # create gwf model
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name,
+                               model_nam_file='{}.nam'.format(name),
+                               save_flows=True)
+
+    # 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',
+                               backtracking_number=0,
+                               inner_maximum=ninner,
+                               inner_hclose=hclose,
+                               rcloserecord=rclose,
+                               linear_acceleration='CG',
+                               scaling_method='NONE',
+                               reordering_method='NONE',
+                               relaxation_factor=relax,
+                               number_orthogonalizations=7)
+    sim.register_ims_package(ims, [gwf.name])
+
+    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
+    flopy.mf6.ModflowGwfic(gwf, strt=strt,
+                           fname='{}.ic'.format(name))
+
+    # node property flow
+    flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=hk)
+
+    # gwf head observations
+    if binaryobs:
+        obs_recarray = {'head0.obs.bsv': obs_data0,
+                        'head1.obs.bsv': obs_data1}
+    else:
+        obs_recarray = {'head0.obs.csv': obs_data0,
+                        'head1.obs.csv': obs_data1}
+
+    o = flopy.mf6.ModflowUtlobs(gwf, pname='head_obs',
+                                fname='{}.obs'.format(name),
+                                digits=10, print_input=True,
+                                continuous=obs_recarray)
+
+    # chd files
+    flopy.mf6.modflow.ModflowGwfchd(gwf, stress_period_data=cd6)
+
+    # output control
+    flopy.mf6.ModflowGwfoc(gwf,
+                           budget_filerecord='{}.cbc'.format(name),
+                           head_filerecord='{}.hds'.format(name),
+                           headprintrecord=[('COLUMNS', 10, 'WIDTH', 15,
+                                             'DIGITS', 6, 'GENERAL')],
+                           saverecord=[('HEAD', 'LAST')],
+                           printrecord=[('HEAD', 'LAST'),
+                                        ('BUDGET', 'LAST')])
+
+    return sim
+
+
+def get_model(idx, dir):
+    ws = dir
+    # build mf6 with ascii observation output
+    sim = build_mf6(idx, ws, binaryobs=False)
+
+    # build mf6 with binary observation output
+    wsc = os.path.join(ws, 'mf6')
+    mc = build_mf6(idx, wsc, binaryobs=True)
+
+    return sim, mc
+
+
+def build_models():
+    for idx, dir in enumerate(exdirs):
+        sim, mc = get_model(idx, dir)
+        sim.write_simulation()
+        mc.write_simulation()
+        hack_binary_obs(idx, dir)
+    return
+
+
+def hack_binary_obs(idx, dir):
+    name = ex[idx]
+    ws = dir
+    wsc = os.path.join(ws, 'mf6')
+    fname = name + '.obs'
+    fpth = os.path.join(wsc, fname)
+    with open(fpth, 'r') as f:
+        lines = f.readlines()
+    with open(fpth, 'w') as f:
+        for line in lines:
+            line = line.rstrip()
+            if 'BEGIN continuous  FILEOUT' in line:
+                line += '  BINARY'
+            f.write('{}\n'.format(line))
+        f.close()
+    return
+
+def eval_obs(sim):
+    print('evaluating observations...')
+
+    # get results from the observation files
+    pth = sim.simpath
+    files = [fn for fn in os.listdir(pth) if '.csv' in fn]
+    for file in files:
+        pth0 = os.path.join(pth, file)
+        pth1 = os.path.join(pth, 'mf6', file.replace('.csv', '.bsv'))
+        d0 = flopy.utils.Mf6Obs(pth0, isBinary=False).get_data()
+        d1 = flopy.utils.Mf6Obs(pth1, isBinary=True).get_data()
+        names0 = d0.dtype.names
+        names1 = d1.dtype.names
+        msg = 'The number of columns ({}) '.format(len(names0)) + \
+              'in {} '.format(pth0) + 'is not equal to ' + \
+              'the number of columns ({}) '.format(len(names1)) + \
+              'in {}.'.format(pth1)
+        assert len(names0) == len(names1), msg
+        msg = 'The number of rows ({}) '.format(d0.shape[0]) + \
+              'in {} '.format(pth0) + 'is not equal to ' + \
+              'the number of rows ({}) '.format(d1.shape[0]) + \
+              'in {}.'.format(pth1)
+        assert d0.shape[0] == d1.shape[0], msg
+        for name in names0:
+            msg = "The values for column '{}' ".format(name) + \
+                  'are not within 1e-5 of each other'
+            assert np.allclose(d0[name], d1[name], rtol=1e-5), 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 dir in exdirs:
+        yield test.run_mf6, Simulation(dir, exfunc=eval_obs)
+
+    return
+
+
+def main():
+    # initialize testing framework
+    test = testing_framework()
+
+    # build the models
+    build_models()
+
+    # run the test models
+    for dir in exdirs:
+        sim = Simulation(dir, exfunc=eval_obs)
+        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

@@ -18,9 +18,9 @@
             "email": "langevin@usgs.gov"
         }, 
         "laborHours": -1, 
-        "version": "6.0.3.30", 
+        "version": "6.0.3.31", 
         "date": {
-            "metadataLastUpdated": "2018-10-29"
+            "metadataLastUpdated": "2018-10-30"
         }, 
         "organization": "U.S. Geological Survey", 
         "permissions": {

doc/ReleaseNotes/ReleaseNotes.tex

@@ -122,7 +122,8 @@ \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 startime 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 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. 
 \end{itemize}
 
 \underline{STRESS PACKAGES}
@@ -141,9 +142,9 @@ \section{History}
 
 \underline{SOLUTION}
 \begin{itemize}
-\item -
-\item -
-\item - 
+\item Modified pseudo-transient continuation (PTC) approach to use PTC for steady-state stress period for models using the Newton-Raphson formulation for problems with and without the storage (STO) package. Previously, PTC was only used with problems that did not include the STO package (this was not the intended behavior of PTC).
+\item Added NO\_PTC option to disable PTC for problems where PTC degrades/prevents model convergence. Option only applies to steady-state stress periods for models using the Newton-Raphson formulation. For many problems, PTC can significantly improve convergence behavior for steady-state simulations, and for this reason it is active by default.  In some cases, however, PTC can worsen the convergence behavior, especially when the initial conditions are similar to the solution.  When the initial conditions are similar to, or exactly the same as, the solution and convergence is slow, then this NO\_PTC option should be used to deactivate PTC.  This NO\_PTC option should also be used in order to compare convergence behavior with other MODFLOW versions, as PTC is only available in MODFLOW 6. 
+\item Small improvements to PTC to reduce the initial PTCDEL value loaded on the diagonal. This reduces the number of iterations required to achieve convergence for steady-state stress periods for most problems.
 \end{itemize}
 
 \item Version mf6.0.3--Aug. 9, 2018

doc/version.tex

@@ -1,3 +1,3 @@
-\newcommand{\modflowversion}{mf6.0.3.30}
-\newcommand{\modflowdate}{October 29, 2018}
+\newcommand{\modflowversion}{mf6.0.3.31}
+\newcommand{\modflowdate}{October 30, 2018}
 \newcommand{\currentmodflowversion}{Version \modflowversion---\modflowdate}

src/Utilities/Observation/Obs3.f90

@@ -846,8 +846,7 @@ subroutine build_headers(this)
     ! -- local
     integer(I4B) :: i, ii, idx, indx, iu, num, nunit
     integer(int32) :: nobs
-    character(len=LENBIGLINE)          :: oldheader, newheader
-    character(len=LENBIGLINE), pointer :: headr => null()
+    character(len=LENOBSNAME), pointer :: headr => null()
     character(len=LENOBSNAME)          :: nam
     character(len=4)                   :: clenobsname
     type(ObserveType),         pointer :: obsrv => null()
@@ -866,13 +865,10 @@ subroutine build_headers(this)
         if (headr == '') then
           headr = 'time'
         endif
-        oldheader = headr
         nam = obsrv%Name
         call ExpandArray(obsOutput%obsnames)
         idx = size(obsOutput%obsnames)
         obsOutput%obsnames(idx) = nam
-        newheader = trim(oldheader) // ',' // trim(nam)
-        headr = newheader
       enddo
     endif
     !
@@ -885,8 +881,13 @@ subroutine build_headers(this)
         ! -- write header to formatted file
         headr => obsOutput%header
         if (headr /= '') then
+          nobs = obsOutput%nobs
           iu = obsOutput%nunit
-          write(iu,'(a)')trim(headr)
+          write(iu, '(a)', advance='NO') 'time'
+          do ii = 1,nobs
+            write(iu, '(a,a)', advance='NO') ',', trim(obsOutput%obsnames(ii))
+          enddo
+          write(iu, '(a)', advance='YES') ''
         endif
       else
         ! -- write header to unformatted file

src/Utilities/Observation/ObsOutput.f90

@@ -25,8 +25,8 @@ module ObsOutputModule
     integer(I4B), public :: nunit = 0
     character(len=500), public :: filename = ''
     character(len=LENOBSNAME), allocatable, dimension(:), public :: obsnames
-    character(len=LENBIGLINE), public :: header = ''
-    character(len=LENBIGLINE), public :: lineout = ''
+    character(len=LENOBSNAME), public :: header = ''
+    character(len=LENOBSNAME), public :: lineout = ''
     logical, public :: FormattedOutput = .true.
   contains
     ! -- Public procedures
@@ -64,8 +64,9 @@ subroutine WriteLineout(this)
     implicit none
     ! -- dummy
     class(ObsOutputType), intent(inout) :: this
-    !
-    write(this%nunit,'(a)')trim(this%lineout)
+    ! -- write a line return to end of observation output line 
+    !    for this totim
+    write(this%nunit,'(a)', advance='YES') ''
     !
     return
   end subroutine WriteLineout