diff --git a/autotest/simulation.py b/autotest/simulation.py index d9eb0b7cdb6..18fb064f30c 100644 --- a/autotest/simulation.py +++ b/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} " diff --git a/autotest/test_z01_testmodels_mf6.py b/autotest/test_z01_testmodels_mf6.py index ad4f6fae38b..de6a5339e94 100644 --- a/autotest/test_z01_testmodels_mf6.py +++ b/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 diff --git a/doc/ReleaseNotes/v6.5.0.tex b/doc/ReleaseNotes/v6.5.0.tex index dde1bd149b9..7abe392aaf6 100644 --- a/doc/ReleaseNotes/v6.5.0.tex +++ b/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} diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90 index b46efce652b..25e4217d56a 100644 --- a/src/Model/GroundWaterFlow/gwf3lak8.f90 +++ b/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 diff --git a/src/Model/GroundWaterFlow/gwf3sfr8.f90 b/src/Model/GroundWaterFlow/gwf3sfr8.f90 index acd52c69f23..21e1c67e65f 100644 --- a/src/Model/GroundWaterFlow/gwf3sfr8.f90 +++ b/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,6 +2206,12 @@ 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 ! @@ -2202,6 +2219,11 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) 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 @@ -2209,7 +2231,15 @@ subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) 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 diff --git a/src/Model/GroundWaterFlow/gwf3uzf8.f90 b/src/Model/GroundWaterFlow/gwf3uzf8.f90 index d3bf9d9022d..5327d5ee779 100644 --- a/src/Model/GroundWaterFlow/gwf3uzf8.f90 +++ b/src/Model/GroundWaterFlow/gwf3uzf8.f90 @@ -1198,9 +1198,12 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) integer(I4B) :: locdrejinfmax integer(I4B) :: locdrchmax integer(I4B) :: locdseepmax + integer(I4B) :: locdqfrommvrmax integer(I4B) :: ntabrows integer(I4B) :: ntabcols integer(I4B) :: n + real(DP) :: q + real(DP) :: q0 real(DP) :: qtolfact real(DP) :: drejinf real(DP) :: drejinfmax @@ -1208,6 +1211,8 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) real(DP) :: drchmax real(DP) :: dseep real(DP) :: dseepmax + real(DP) :: dqfrommvr + real(DP) :: dqfrommvrmax ! format ! -------------------------------------------------------------------------- ! @@ -1217,9 +1222,11 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) locdrejinfmax = 0 locdrchmax = 0 locdseepmax = 0 + locdqfrommvrmax = 0 drejinfmax = DZERO drchmax = DZERO dseepmax = DZERO + dqfrommvrmax = DZERO ! ! -- if not saving package convergence data on check convergence if ! the model is considered converged @@ -1238,6 +1245,9 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) if (this%iseepflag == 1) 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, '') @@ -1270,6 +1280,12 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) tag = 'dseepmax_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 ! @@ -1292,6 +1308,14 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) dseep = qtolfact * (this%gwd0(n) - this%gwd(n)) 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 drejinfmax = drejinf @@ -1300,6 +1324,8 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) locdrchmax = n dseepmax = dseep locdseepmax = n + dqfrommvrmax = dqfrommvr + locdqfrommvrmax = n else if (ABS(drejinf) > abs(drejinfmax)) then drejinfmax = drejinf @@ -1313,6 +1339,10 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) dseepmax = dseep locdseepmax = n end if + if (ABS(dqfrommvr) > abs(dqfrommvrmax)) then + dqfrommvrmax = dqfrommvr + locdqfrommvrmax = n + end if end if end do final_check ! @@ -1337,6 +1367,14 @@ subroutine uzf_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 @@ -1355,6 +1393,10 @@ subroutine uzf_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) call this%pakcsvtab%add_term(dseepmax) call this%pakcsvtab%add_term(locdseepmax) 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 diff --git a/src/Model/ModelUtilities/PackageMover.f90 b/src/Model/ModelUtilities/PackageMover.f90 index b0c43187838..98dc2a9c625 100644 --- a/src/Model/ModelUtilities/PackageMover.f90 +++ b/src/Model/ModelUtilities/PackageMover.f90 @@ -17,10 +17,11 @@ module PackageMoverModule integer(I4B), pointer :: nproviders integer(I4B), pointer :: nreceivers integer(I4B), dimension(:), pointer, contiguous :: iprmap => null() !< map between id1 and feature (needed for lake to map from outlet to lake number) - real(DP), dimension(:), pointer, contiguous :: qtformvr => null() - real(DP), dimension(:), pointer, contiguous :: qformvr => null() - real(DP), dimension(:), pointer, contiguous :: qtomvr => null() - real(DP), dimension(:), pointer, contiguous :: qfrommvr => null() + real(DP), dimension(:), pointer, contiguous :: qtformvr => null() !< total flow rate available for mover + real(DP), dimension(:), pointer, contiguous :: qformvr => null() !< currently available consumed water (changes during fc) + real(DP), dimension(:), pointer, contiguous :: qtomvr => null() !< actual amount of water sent to mover + real(DP), dimension(:), pointer, contiguous :: qfrommvr => null() !< actual amount of water received from mover + real(DP), dimension(:), pointer, contiguous :: qfrommvr0 => null() !< qfrommvr from previous iteration contains procedure :: ar @@ -31,6 +32,7 @@ module PackageMoverModule procedure :: allocate_scalars procedure :: allocate_arrays procedure :: get_qfrommvr + procedure :: get_qfrommvr0 procedure :: get_qtomvr procedure :: accumulate_qformvr @@ -51,6 +53,7 @@ subroutine set_packagemover_pointer(packagemover, memPath) call mem_setptr(packagemover%qformvr, 'QFORMVR', packagemover%memoryPath) call mem_setptr(packagemover%qtomvr, 'QTOMVR', packagemover%memoryPath) call mem_setptr(packagemover%qfrommvr, 'QFROMMVR', packagemover%memoryPath) + call mem_setptr(packagemover%qfrommvr0, 'QFROMMVR0', packagemover%memoryPath) end subroutine set_packagemover_pointer subroutine nulllify_packagemover_pointer(packagemover) @@ -63,6 +66,7 @@ subroutine nulllify_packagemover_pointer(packagemover) packagemover%qformvr => null() packagemover%qtomvr => null() packagemover%qfrommvr => null() + packagemover%qfrommvr0 => null() end subroutine nulllify_packagemover_pointer subroutine ar(this, nproviders, nreceivers, memoryPath) @@ -102,6 +106,7 @@ subroutine cf(this) ! ! -- set frommvr and qtomvr to zero do i = 1, this%nreceivers + this%qfrommvr0(i) = this%qfrommvr(i) this%qfrommvr(i) = DZERO end do do i = 1, this%nproviders @@ -135,6 +140,7 @@ subroutine da(this) call mem_deallocate(this%qformvr) call mem_deallocate(this%qtomvr) call mem_deallocate(this%qfrommvr) + call mem_deallocate(this%qfrommvr0) ! ! -- scalars call mem_deallocate(this%nproviders) @@ -169,6 +175,8 @@ subroutine allocate_arrays(this) call mem_allocate(this%qformvr, this%nproviders, 'QFORMVR', this%memoryPath) call mem_allocate(this%qtomvr, this%nproviders, 'QTOMVR', this%memoryPath) call mem_allocate(this%qfrommvr, this%nreceivers, 'QFROMMVR', this%memoryPath) + call mem_allocate(this%qfrommvr0, this%nreceivers, 'QFROMMVR0', & + this%memoryPath) ! ! -- initialize do i = 1, this%nproviders @@ -179,6 +187,7 @@ subroutine allocate_arrays(this) end do do i = 1, this%nreceivers this%qfrommvr(i) = DZERO + this%qfrommvr0(i) = DZERO end do ! ! -- return @@ -192,6 +201,13 @@ function get_qfrommvr(this, ireceiver) result(qfrommvr) qfrommvr = this%qfrommvr(ireceiver) end function get_qfrommvr + function get_qfrommvr0(this, ireceiver) result(qfrommvr0) + class(PackageMoverType) :: this + real(DP) :: qfrommvr0 + integer, intent(in) :: ireceiver + qfrommvr0 = this%qfrommvr0(ireceiver) + end function get_qfrommvr0 + function get_qtomvr(this, iprovider) result(qtomvr) class(PackageMoverType) :: this real(DP) :: qtomvr