fix(src_cc): Add SFR convergence check (#85)

Add OUTER_RCLOSEBND variable to IMS that is used when performing final
convergence checks on model packages that solve a separate equation
not solved by the IMS linear solver. Add stage and residual convergence
checks to the SFR package to make sure that stage and upstream flow
changes between successive outer iterations are less than OUTER_HCLOSE
and OUTER_RCLOSEBND, respectively.

Modify the final convergence check for the LAK package to use
OUTER_HCLOSE when evaluating lake stage changes between successive
outer iterations. Modify the final convergence check for the UZF package
to use OUTER_RCLOSEBND when evaluating rejected infiltration,
groundwater recharge, and groundwater seepage changes between successive
outer iterations.

README.md

@@ -7,7 +7,7 @@
 
 ## Automated Testing Status on Travis-CI
 
-### Version 6.0.3 develop — build 56
+### Version 6.0.3 develop — build 42
 [![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 February 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, 08 February 2019, https://doi.org/10.5066/F76Q1VQV](https://doi.org/10.5066/F76Q1VQV)
 
 ## Instructions for building definition files for new packages
 

code.json

@@ -18,9 +18,9 @@
             "email": "langevin@usgs.gov"
         }, 
         "laborHours": -1, 
-        "version": "6.0.3.56", 
+        "version": "6.0.3.42", 
         "date": {
-            "metadataLastUpdated": "2019-02-07"
+            "metadataLastUpdated": "2019-02-08"
         }, 
         "organization": "U.S. Geological Survey", 
         "permissions": {

doc/ReleaseNotes/ReleaseNotes.tex

@@ -137,7 +137,10 @@ \section{History}
 \underline{ADVANCED STRESS PACKAGES}
 \begin{itemize}
 \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 Added stage and residual convergence checks to the SFR package to make sure that stage and upstream flow changes between successive outer iterations are less than OUTER\_HCLOSE and OUTER\_RCLOSEBND, respectively. This addition is expected to be useful for steady-state simulations with complicated networks and simple reaches.
+\item Modified the final convergence check for the LAK package to use OUTER\_HCLOSE when evaluating lake stage changes between successive outer iterations.
+\item Modified the final convergence check for the UZF package to use OUTER\_RCLOSEBND when evaluating rejected infiltration, groundwater recharge, and groundwater seepage changes between successive outer iterations.
+\item - 
 \item - 
 \end{itemize}
 
@@ -146,6 +149,8 @@ \section{History}
 \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.
+\item Added OUTER\_RCLOSEBND variable that is used when performing final convergence checks on model packages that solve a separate equation not solved by the IMS linear solver. This value represents the maximum allowable residual at any single model package element between successive outer iterations. An example of a model package that would use OUTER\_RCLOSEBND to evaluate convergence is the SFR package which solves a continuity equation for each reach.
+\item - 
 \end{itemize}
 
 \item Version mf6.0.3--Aug. 9, 2018

doc/mf6io/ims_table.tex

@@ -4,6 +4,7 @@
 Nonlinear Variable & default/\texttt{simple} & \texttt{moderate} & \texttt{complex} \\
 \hline
 \texttt{OUTER\_HCLOSE} & 0.001 & 0.01 & 0.1 \\
+\texttt{OUTER\_RCLOSEBND} & 0.1 & 0.1 & 0.1 \\
 \texttt{OUTER\_MAXIMUM} & 25 & 50 & 100 \\
 \texttt{UNDER\_RELAXATION} & \texttt{NONE} & \texttt{DBD} & \texttt{DBD} \\
 \texttt{UNDER\_RELAXATION\_THETA} & 0.0 & 0.9 & 0.8 \\

doc/mf6io/mf6ivar/dfn/sln-ims.dfn

@@ -79,6 +79,14 @@ optional false
 longname head change criterion
 description real value defining the head change criterion for convergence of the outer (nonlinear) iterations, in units of length. When the maximum absolute value of the head change at all nodes during an iteration is less than or equal to OUTER\_HCLOSE, iteration stops. Commonly, OUTER\_HCLOSE equals 0.01.
 
+block nonlinear
+name outer_rclosebnd
+type double precision
+reader urword
+optional true
+longname boundary package flow residual tolerance
+description real value defining the residual tolerance for convergence of model packages that solve a separate equation not solved by the IMS linear solver. This value represents the maximum allowable residual between successive outer iterations at any single model package element. An example of a model package that would use OUTER\_RCLOSEBND to evaluate convergence is the SFR package which solves a continuity equation for each reach.   
+
 block nonlinear
 name outer_maximum
 type integer

doc/mf6io/mf6ivar/md/mf6ivar.md

@@ -54,6 +54,7 @@
 | SLN | IMS | OPTIONS | CSVFILE | STRING | name of the ascii comma separated values output file to write solver convergence information. If PRINT\_OPTION is NONE or SUMMARY, comma separated values output includes maximum head change convergence information at the end of each outer iteration for each time step. If PRINT\_OPTION is ALL, comma separated values output includes maximum head change and maximum residual convergence information for the solution and each model (if the solution includes more than one model) and linear acceleration information for each inner iteration. |
 | SLN | IMS | OPTIONS | NO_PTC | KEYWORD | is a flag that is used to disable pseudo-transient continuation (PTC). 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. |
 | SLN | IMS | NONLINEAR | OUTER_HCLOSE | DOUBLE PRECISION | real value defining the head change criterion for convergence of the outer (nonlinear) iterations, in units of length. When the maximum absolute value of the head change at all nodes during an iteration is less than or equal to OUTER\_HCLOSE, iteration stops. Commonly, OUTER\_HCLOSE equals 0.01. |
+| SLN | IMS | NONLINEAR | OUTER_RCLOSEBND | DOUBLE PRECISION | real value defining the residual tolerance for convergence of model packages that solve a separate equation not solved by the IMS linear solver. This value represents the maximum allowable residual between successive outer iterations at any single model package element. An example of a model package that would use OUTER\_RCLOSEBND to evaluate convergence is the SFR package which solves a continuity equation for each reach. |
 | SLN | IMS | NONLINEAR | OUTER_MAXIMUM | INTEGER | integer value defining the maximum number of outer (nonlinear) iterations -- that is, calls to the solution routine. For a linear problem OUTER\_MAXIMUM should be 1. |
 | SLN | IMS | NONLINEAR | UNDER_RELAXATION | STRING | is an optional keyword that defines the nonlinear under-relaxation schemes used. By default under-relaxation is not used.  NONE - under-relaxation is not used. SIMPLE - Simple under-relaxation scheme with a fixed relaxation factor is used.  COOLEY - Cooley under-relaxation scheme is used.  DBD - delta-bar-delta under-relaxation is used.  Note that the under-relaxation schemes are used in conjunction with problems that use the Newton-Raphson formulation, however, experience has indicated that the Cooley under-relaxation and damping work well also for the Picard scheme with the wet/dry options of MODFLOW 6. |
 | SLN | IMS | NONLINEAR | UNDER_RELAXATION_THETA | DOUBLE PRECISION | real value defining the reduction factor for the learning rate (under-relaxation term) of the delta-bar-delta algorithm. The value of UNDER\_RELAXATION\_THETA is between zero and one. If the change in the variable (head) is of opposite sign to that of the previous iteration, the under-relaxation term is reduced by a factor of UNDER\_RELAXATION\_THETA. The value usually ranges from 0.3 to 0.9; a value of 0.7 works well for most problems. UNDER\_RELAXATION\_THETA only needs to be specified if UNDER\_RELAXATION is DBD. |

doc/mf6io/mf6ivar/tex/sln-ims-desc.tex

@@ -21,6 +21,8 @@
 \begin{description}
 \item \texttt{outer\_hclose}---real value defining the head change criterion for convergence of the outer (nonlinear) iterations, in units of length. When the maximum absolute value of the head change at all nodes during an iteration is less than or equal to OUTER\_HCLOSE, iteration stops. Commonly, OUTER\_HCLOSE equals 0.01.
 
+\item \texttt{outer\_rclosebnd}---real value defining the residual tolerance for convergence of model packages that solve a separate equation not solved by the IMS linear solver. This value represents the maximum allowable residual between successive outer iterations at any single model package element. An example of a model package that would use OUTER\_RCLOSEBND to evaluate convergence is the SFR package which solves a continuity equation for each reach.
+
 \item \texttt{outer\_maximum}---integer value defining the maximum number of outer (nonlinear) iterations -- that is, calls to the solution routine. For a linear problem OUTER\_MAXIMUM should be 1.
 
 \item \texttt{under\_relaxation}---is an optional keyword that defines the nonlinear under-relaxation schemes used. By default under-relaxation is not used.  NONE - under-relaxation is not used. SIMPLE - Simple under-relaxation scheme with a fixed relaxation factor is used.  COOLEY - Cooley under-relaxation scheme is used.  DBD - delta-bar-delta under-relaxation is used.  Note that the under-relaxation schemes are used in conjunction with problems that use the Newton-Raphson formulation, however, experience has indicated that the Cooley under-relaxation and damping work well also for the Picard scheme with the wet/dry options of MODFLOW 6.

doc/mf6io/mf6ivar/tex/sln-ims-nonlinear.dat

@@ -1,5 +1,6 @@
 BEGIN NONLINEAR
   OUTER_HCLOSE <outer_hclose>
+  [OUTER_RCLOSEBND <outer_rclosebnd>]
   OUTER_MAXIMUM <outer_maximum>
   [UNDER_RELAXATION <under_relaxation>]
   [UNDER_RELAXATION_THETA <under_relaxation_theta>]

doc/version.tex

@@ -1,3 +1,3 @@
-\newcommand{\modflowversion}{mf6.0.3.56}
-\newcommand{\modflowdate}{February 07, 2019}
+\newcommand{\modflowversion}{mf6.0.3.42}
+\newcommand{\modflowdate}{February 08, 2019}
 \newcommand{\currentmodflowversion}{Version \modflowversion---\modflowdate}

src/Model/GroundWaterFlow/gwf3.f90

@@ -652,7 +652,7 @@ subroutine gwf_fc(this, kiter, amatsln, njasln, inwtflag)
     return
   end subroutine gwf_fc
 
-  subroutine gwf_cc(this, kiter, iend, icnvg)
+  subroutine gwf_cc(this, kiter, iend, icnvg, hclose, rclose)
 ! ******************************************************************************
 ! gwf_cc -- GroundWater Flow Model Final Convergence Check for Boundary Packages
 ! Subroutine: (1) calls package cc routines
@@ -665,6 +665,8 @@ subroutine gwf_cc(this, kiter, iend, icnvg)
     integer(I4B),intent(in) :: kiter
     integer(I4B),intent(in) :: iend
     integer(I4B),intent(inout) :: icnvg
+    real(DP), intent(in) :: hclose
+    real(DP), intent(in) :: rclose
     ! -- local
     class(BndType), pointer :: packobj
     integer(I4B) :: ip
@@ -677,7 +679,7 @@ subroutine gwf_cc(this, kiter, iend, icnvg)
     ! -- Call package cc routines
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
-      call packobj%bnd_cc(iend, icnvg)
+      call packobj%bnd_cc(iend, icnvg, hclose, rclose)
     enddo
     !
     ! -- return

src/Model/GroundWaterFlow/gwf3lak8.f90

@@ -2350,8 +2350,8 @@ end subroutine lak_estimate_conn_exchange
 
   subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
 ! ******************************************************************************
-! lak_calculate_storagechange -- Calculate the inflow terms to a lake at a
-!                         provided stage.
+! lak_calculate_storagechange -- Calculate the storage change in a lake based on
+!                         provided stages and a passed delt.
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
@@ -3786,22 +3786,27 @@ subroutine lak_fn(this, rhs, ia, idxglo, amatsln)
     return
   end subroutine lak_fn
 
-  subroutine lak_cc(this, iend, icnvg)
+  subroutine lak_cc(this, iend, icnvg, hclose, rclose)
 ! **************************************************************************
 ! lak_cc -- Final convergence check for package
 ! **************************************************************************
 !
 !    SPECIFICATIONS:
 ! --------------------------------------------------------------------------
+    use TdisModule, only: delt
     ! -- dummy
     class(LakType), intent(inout) :: this
     integer(I4B), intent(in) :: iend
     integer(I4B), intent(inout) :: icnvg
+    real(DP), intent(in) :: hclose
+    real(DP), intent(in) :: rclose
     ! -- local
     integer(I4B) :: n
     integer(I4B) :: ifirst
     real(DP) :: dh
+    real(DP) :: residb0
     real(DP) :: residb
+    real(DP) :: dr
     real(DP) :: inf
     real(DP) :: outf
     real(DP) :: avgf
@@ -3821,8 +3826,10 @@ subroutine lak_cc(this, iend, icnvg)
     if (this%iconvchk /= 0) then
       final_check: do n = 1, this%nlakes
         if (this%iboundpak(n) < 1) cycle
-        dh = ABS(this%s0(n) - this%xnewpak(n))
+        dh = this%s0(n) - this%xnewpak(n)
+        call this%lak_calculate_residual(n, this%s0(n), residb0)
         call this%lak_calculate_residual(n, this%xnewpak(n), residb)
+        dr = residb0 - residb
         call this%lak_calculate_available(n, this%xnewpak(n), inf, &
                                           ra, ro, qinf, ex)
         outf = inf - residb
@@ -3834,18 +3841,19 @@ subroutine lak_cc(this, iend, icnvg)
           end if
         end if
         !write(*,'(1x,i4,6(1x,g10.4))') n, this%s0(n), this%xnewpak(n), residb, outf, inf, pd
-        if (dh > this%delh .or. ABS(pd) > this%pdmax) then
+        if (ABS(dh) > hclose .or. ABS(pd) > this%pdmax) then
           icnvg = 0
           ! write convergence check information if this is the last outer iteration
           if (iend == 1) then
             if (ifirst == 1) then
               ifirst = 0
               write(*,2030) this%name
               write(this%iout, 2000) '      LAKE',                                 &
-                '        MAX. DH', '    DH CRITERIA',                              &
-                '      PCT DIFF.', 'PCT DIFF. CRIT.'
+                '             DH', '    DH CRITERIA',                              &
+                '        PCTDIFF', ' PCTDIFF CRITER'
             end if
-            write(this%iout,2010) n, dh, this%delh, pd, this%pdmax
+            !write(this%iout,2010) n, dh, this%delh, pd, this%pdmax
+            write(this%iout,2010) n, dh, hclose, pd, this%pdmax
           else
             exit final_check
           end if
@@ -6020,29 +6028,6 @@ subroutine lak_calculate_residual(this, n, hlak, resid, headp)
     ! -- calculate the available water
     call this%lak_calculate_available(n, hlak, avail, &
                                       ra, ro, qinf, ex, hp)
-    !!
-    !! -- calculate the aquifer sources to the lake
-    !do j = this%idxlakeconn(n), this%idxlakeconn(n+1)-1
-    !  igwfnode = this%cellid(j)
-    !  head = this%xnew(igwfnode) + hp
-    !  call this%lak_estimate_conn_exchange(1, n, j, idry, hlak, head, qlakgw, clak, avail)
-    !end do
-    !!
-    !! -- add rainfall
-    !call this%lak_calculate_rainfall(n, hlak, ra)
-    !avail = avail + ra
-    !!
-    !! -- calculate runoff
-    !call this%lak_calculate_runoff(n, ro)
-    !avail = avail + ro
-    !!
-    !! -- calculate inflow
-    !call this%lak_calculate_inflow(n, qinf)
-    !avail = avail + qinf
-    !!
-    !! -- calculate external flow terms
-    !call this%lak_calculate_external(n, ex)
-    !avail = avail + ex
     !
     ! -- calculate groundwater seepage
     do j = this%idxlakeconn(n), this%idxlakeconn(n+1)-1