feat(uzf): add qfrommvr convergence check to UZF, LAK, SFR (#1119)

autotest/simulation.py

@@ -771,7 +771,7 @@ def _compare_budgets(self, msgall, extensions="cbc"):
                             f"{os.path.basename(fpth0)} - "
                             + f"{key:16s} "
                             + f"difference ({diffmax:10.4g}) "
-                            + f"> {self.pdtol:10.4g} "
+                            + f"> {vmin_tol:10.4g} "
                             + f"at {indices.size} nodes "
                             + f" [first location ({indices[0] + 1})] "
                             + f"at time {t} "

autotest/test_z01_testmodels_mf6.py

@@ -189,8 +189,8 @@ def set_make_comparison(test):
         "test051_uzfp3_lakmvr_v2": ("6.2.1",),
         "test051_uzfp3_wellakmvr_v2": ("6.2.1",),
         "test045_lake4ss": ("6.2.2",),
-        "test056_mt3dms_usgs_gwtex_dev": ("6.2.2",),
-        "test056_mt3dms_usgs_gwtex_IR_dev": ("6.2.2",),
+        "test056_mt3dms_usgs_gwtex_dev": ("6.4.1",),
+        "test056_mt3dms_usgs_gwtex_IR_dev": ("6.4.1",),
     }
     make_comparison = True
     if test in compare_tests.keys():
@@ -200,6 +200,7 @@ def set_make_comparison(test):
         print(f"MODFLOW regression version='{version}'")
         if version in compare_tests[test]:
             make_comparison = False
+            print(f"Make comparison has been set to False.")
     return make_comparison
 
 

doc/ReleaseNotes/v6.5.0.tex

@@ -39,12 +39,12 @@
 	%	\item xxx
 	%\end{itemize}
 
-	%\underline{ADVANCED STRESS PACKAGES}
-	%\begin{itemize}
-	%	\item xxx
+	\underline{ADVANCED STRESS PACKAGES}
+	\begin{itemize}
+		\item Added additional convergence checks to the Streamflow Routing (SFR), Lake (LAK) and Unsaturated Zone Flow (UZF) Packages to ensure that flows from the Water Mover (MVR) Package meet solver tolerance.  Mover flows are converted into depths using the time step length and area, and the depths are compared to the Iterative Model Solution (IMS) DVCLOSE input parameter.  If a depth is greater than DVCLOSE, then the iteration is marked as not converged.  The maximum depth change between iterations and the advanced package feature number is written for each outer iteration to a comma-separated value file, provided the PACKAGE\_CONVERGENCE option is specified in the options block.
 	%	\item xxx
 	%	\item xxx
-	%\end{itemize}
+	\end{itemize}
 
 	%\underline{SOLUTION}
 	%\begin{itemize}

src/Model/GroundWaterFlow/gwf3lak8.f90

@@ -4024,9 +4024,13 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
     integer(I4B) :: locdhmax
     integer(I4B) :: locdgwfmax
     integer(I4B) :: locdqoutmax
+    integer(I4B) :: locdqfrommvrmax
     integer(I4B) :: ntabrows
     integer(I4B) :: ntabcols
     integer(I4B) :: n
+    real(DP) :: q
+    real(DP) :: q0
+    real(DP) :: qtolfact
     real(DP) :: area
     real(DP) :: gwf0
     real(DP) :: gwf
@@ -4045,6 +4049,8 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
     real(DP) :: dhmax
     real(DP) :: dgwfmax
     real(DP) :: dqoutmax
+    real(DP) :: dqfrommvr
+    real(DP) :: dqfrommvrmax
     ! format
 ! --------------------------------------------------------------------------
     !
@@ -4054,9 +4060,11 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
     locdhmax = 0
     locdgwfmax = 0
     locdqoutmax = 0
+    locdqfrommvrmax = 0
     dhmax = DZERO
     dgwfmax = DZERO
     dqoutmax = DZERO
+    dqfrommvrmax = DZERO
     !
     ! -- if not saving package convergence data on check convergence if
     !    the model is considered converged
@@ -4077,6 +4085,9 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         if (this%noutlets > 0) then
           ntabcols = ntabcols + 2
         end if
+        if (this%imover == 1) then
+          ntabcols = ntabcols + 2
+        end if
         !
         ! -- setup table
         call table_cr(this%pakcsvtab, this%packName, '')
@@ -4109,6 +4120,12 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
           tag = 'dqoutmax_loc'
           call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT)
         end if
+        if (this%imover == 1) then
+          tag = 'dqfrommvrmax'
+          call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT)
+          tag = 'dqfrommvrmax_loc'
+          call this%pakcsvtab%initialize_column(tag, 16, alignment=TABLEFT)
+        end if
       end if
     end if
     !
@@ -4127,12 +4144,15 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         ! -- calculate surface area
         call this%lak_calculate_sarea(n, hlak, area)
         !
+        ! -- set the Q to length factor
+        qtolfact = delt / area
+        !
         ! -- change in gwf exchange
         dgwf = DZERO
         if (area > DZERO) then
           gwf0 = this%qgwf0(n)
           call this%lak_calculate_exchange(n, hlak, gwf)
-          dgwf = (gwf0 - gwf) * delt / area
+          dgwf = (gwf0 - gwf) * qtolfact
         end if
         !
         ! -- change in outflows
@@ -4143,10 +4163,18 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
             call this%lak_calculate_outlet_outflow(n, hlak0, inf, qout0)
             call this%lak_calculate_available(n, hlak, inf, ra, ro, qinf, ex)
             call this%lak_calculate_outlet_outflow(n, hlak, inf, qout)
-            dqout = (qout0 - qout) * delt / area
+            dqout = (qout0 - qout) * qtolfact
           end if
         end if
         !
+        ! -- q from mvr
+        dqfrommvr = DZERO
+        if (this%imover == 1) then
+          q = this%pakmvrobj%get_qfrommvr(n)
+          q0 = this%pakmvrobj%get_qfrommvr0(n)
+          dqfrommvr = qtolfact * (q0 - q)
+        end if
+        !
         ! -- evaluate magnitude of differences
         if (n == 1) then
           locdhmax = n
@@ -4155,6 +4183,8 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
           dgwfmax = dgwf
           locdqoutmax = n
           dqoutmax = dqout
+          dqfrommvrmax = dqfrommvr
+          locdqfrommvrmax = n
         else
           if (abs(dh) > abs(dhmax)) then
             locdhmax = n
@@ -4168,6 +4198,10 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
             locdqoutmax = n
             dqoutmax = dqout
           end if
+          if (ABS(dqfrommvr) > abs(dqfrommvrmax)) then
+            dqfrommvrmax = dqfrommvr
+            locdqfrommvrmax = n
+          end if
         end if
       end do final_check
       !
@@ -4195,6 +4229,14 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
           cpak = trim(cloc)
         end if
       end if
+      if (this%imover == 1) then
+        if (ABS(dqfrommvrmax) > abs(dpak)) then
+          ipak = locdqfrommvrmax
+          dpak = dqfrommvrmax
+          write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
+          cpak = trim(cloc)
+        end if
+      end if
       !
       ! -- write convergence data to package csv
       if (this%ipakcsv /= 0) then
@@ -4213,6 +4255,10 @@ subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
           call this%pakcsvtab%add_term(dqoutmax)
           call this%pakcsvtab%add_term(locdqoutmax)
         end if
+        if (this%imover == 1) then
+          call this%pakcsvtab%add_term(dqfrommvrmax)
+          call this%pakcsvtab%add_term(locdqfrommvrmax)
+        end if
         !
         ! -- finalize the package csv
         if (iend == 1) then

src/Model/GroundWaterFlow/gwf3sfr8.f90

@@ -2136,13 +2136,19 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
     integer(I4B) :: ipakfail
     integer(I4B) :: locdhmax
     integer(I4B) :: locrmax
+    integer(I4B) :: locdqfrommvrmax
     integer(I4B) :: ntabrows
     integer(I4B) :: ntabcols
     integer(I4B) :: n
+    real(DP) :: q
+    real(DP) :: q0
+    real(DP) :: qtolfact
     real(DP) :: dh
     real(DP) :: r
     real(DP) :: dhmax
     real(DP) :: rmax
+    real(DP) :: dqfrommvr
+    real(DP) :: dqfrommvrmax
     !
     ! -- initialize local variables
     icheck = this%iconvchk
@@ -2152,6 +2158,8 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
     r = DZERO
     dhmax = DZERO
     rmax = DZERO
+    locdqfrommvrmax = 0
+    dqfrommvrmax = DZERO
     !
     ! -- if not saving package convergence data on check convergence if
     !    the model is considered converged
@@ -2169,6 +2177,9 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         ! -- determine the number of columns and rows
         ntabrows = 1
         ntabcols = 9
+        if (this%imover == 1) then
+          ntabcols = ntabcols + 2
+        end if
         !
         ! -- setup table
         call table_cr(this%pakcsvtab, this%packName, '')
@@ -2195,21 +2206,40 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT)
         tag = 'dinflowmax_loc'
         call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT)
+        if (this%imover == 1) then
+          tag = 'dqfrommvrmax'
+          call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT)
+          tag = 'dqfrommvrmax_loc'
+          call this%pakcsvtab%initialize_column(tag, 16, alignment=TABLEFT)
+        end if
       end if
     end if
     !
     ! -- perform package convergence check
     if (icheck /= 0) then
       final_check: do n = 1, this%maxbound
         if (this%iboundpak(n) == 0) cycle
+        !
+        ! -- set the Q to length factor
+        qtolfact = delt / this%calc_surface_area(n)
+        !
+        ! -- calculate stage change
         dh = this%stage0(n) - this%stage(n)
         !
         ! -- evaluate flow difference if the time step is transient
         if (this%gwfiss == 0) then
           r = this%usflow0(n) - this%usflow(n)
           !
           ! -- normalize flow difference and convert to a depth
-          r = r * delt / this%calc_surface_area(n)
+          r = r * qtolfact
+        end if
+        !
+        ! -- q from mvr
+        dqfrommvr = DZERO
+        if (this%imover == 1) then
+          q = this%pakmvrobj%get_qfrommvr(n)
+          q0 = this%pakmvrobj%get_qfrommvr0(n)
+          dqfrommvr = qtolfact * (q0 - q)
         end if
         !
         ! -- evaluate magnitude of differences
@@ -2218,6 +2248,8 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
           dhmax = dh
           locrmax = n
           rmax = r
+          dqfrommvrmax = dqfrommvr
+          locdqfrommvrmax = n
         else
           if (abs(dh) > abs(dhmax)) then
             locdhmax = n
@@ -2227,6 +2259,10 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
             locrmax = n
             rmax = r
           end if
+          if (ABS(dqfrommvr) > abs(dqfrommvrmax)) then
+            dqfrommvrmax = dqfrommvr
+            locdqfrommvrmax = n
+          end if
         end if
       end do final_check
       !
@@ -2243,6 +2279,14 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         write (cloc, "(a,'-',a)") trim(this%packName), 'inflow'
         cpak = trim(cloc)
       end if
+      if (this%imover == 1) then
+        if (ABS(dqfrommvrmax) > abs(dpak)) then
+          ipak = locdqfrommvrmax
+          dpak = dqfrommvrmax
+          write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
+          cpak = trim(cloc)
+        end if
+      end if
       !
       ! -- write convergence data to package csv
       if (this%ipakcsv /= 0) then
@@ -2257,6 +2301,10 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
         call this%pakcsvtab%add_term(locdhmax)
         call this%pakcsvtab%add_term(rmax)
         call this%pakcsvtab%add_term(locrmax)
+        if (this%imover == 1) then
+          call this%pakcsvtab%add_term(dqfrommvrmax)
+          call this%pakcsvtab%add_term(locdqfrommvrmax)
+        end if
         !
         ! -- finalize the package csv
         if (iend == 1) then