fix(ats-budget): fix cumulative budget for failed/retried ATS timeste…

…ps (#1714)

* fix(ats-budget): cumulative budget incorrect with failed/retried ATS timesteps

* add and require a finalize_step method for the budget object
* cumulative terms are calculated only when the finalize_step method is called
* required adding calls to budget%finalize_step() before any budget_ot() calls were made

* update release notes

autotest/test_gwf_ats_lak01.py

@@ -5,7 +5,7 @@
 """
 
 import os
-
+import pathlib as pl
 import flopy
 import numpy as np
 import pytest
@@ -212,6 +212,7 @@ def build_models(idx, test):
     oc = flopy.mf6.ModflowGwfoc(
         gwf,
         budget_filerecord=f"{gwfname}.cbc",
+        budgetcsv_filerecord=f"{gwfname}.bud.csv",
         head_filerecord=f"{gwfname}.hds",
         headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
         saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
@@ -290,6 +291,67 @@ def get_kij_from_node(node, nrow, ncol):
     return k, i, j
 
 
+def budcsv_to_cumulative(fpth):
+    budcsv = np.genfromtxt(fpth, names=True, delimiter=",", deletechars="")
+    nrow = budcsv.shape[0]
+    budcsv_cumulative = budcsv.copy()
+    budcsv_cumulative.resize(nrow + 1)
+    budcsv_cumulative["time"][0] = 0.
+    budcsv_cumulative["time"][1:] = budcsv["time"][:]
+    for name in budcsv.dtype.names[1:]:
+        budcsv_cumulative[name][0] = 0.
+        for i in range(nrow):
+            dt = budcsv_cumulative["time"][i + 1] - budcsv_cumulative["time"][i]
+            budcsv_cumulative[name][i + 1] = budcsv_cumulative[name][i] + budcsv[name][i] * dt
+    return budcsv_cumulative
+
+
+def listfile_to_cumulative(listfile):
+    import flopy
+    mflist = flopy.utils.Mf6ListBudget(listfile)
+    return mflist.get_cumulative()
+
+
+def compare_listbudget_and_budgetcsv(listfile, budcsvfile, verbose, check, atol):
+    """Read a budgetcsv file, convert it to a cumulative budget
+    and then compare it with the cumulative budget in a list file"""
+
+    if verbose:
+        print(f"Comparing {listfile} with {budcsvfile}")
+
+    # get a cumulative budget from the budcsv file
+    budcsvcum = budcsv_to_cumulative(budcsvfile)
+
+    # get the cumulative budget from the list file
+    budlstcum = listfile_to_cumulative(listfile)
+
+    # if print budget is not active for every time step, then the list file
+    # budget may not be complete and comparable to budcsvfile
+    assert budcsvcum.shape[0] - 1 == budlstcum.shape[0], "File sizes are different."
+
+    allclose_list = []
+    for name1 in budlstcum.dtype.names[3:]:
+        nl = name1.split("_")
+        if len(nl) > 1:            
+            for name2 in budcsvcum.dtype.names:
+                if nl[0] in name2 and nl[1] in name2:
+                    # print(f"Found match: {name1} and {name2}")
+                    diff = budcsvcum[name2][1:] - budlstcum[name1]
+                    mindiff = diff.min()
+                    maxdiff = diff.max()
+                    allclose = np.allclose(budcsvcum[name2][1:], budlstcum[name1], atol=atol)
+                    msg = f"{name2} is same: {allclose}.  Min diff: {mindiff} Max diff {maxdiff}"
+                    if verbose:
+                        print(msg)
+                    allclose_list.append((allclose, name1, mindiff, maxdiff, msg))
+
+    if check:
+        for rec in allclose_list:
+            assert rec[0], rec[-1]
+
+    return allclose_list
+
+
 def check_output(idx, test):
     # calculate volume of water and make sure it is conserved
     fname = test.name + ".lak.bin"
@@ -407,6 +469,13 @@ def check_output(idx, test):
     # make_plot(sim, times, head, stage)
     # make_plot_xsect(sim, head, stage)
 
+    listfile = pl.Path(test.workspace) / f"{test.name}.lst"
+    budcsvfile = pl.Path(test.workspace) / f"{test.name}.bud.csv"
+    verbose = True
+    check = True
+    atol = 0.001
+    compare_listbudget_and_budgetcsv(listfile, budcsvfile, verbose, check, atol)
+
 
 @pytest.mark.slow
 @pytest.mark.parametrize("idx, name", enumerate(cases))

doc/ReleaseNotes/develop.tex

@@ -18,12 +18,12 @@
 	%\end{itemize}
 
 	\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\
-	%\underline{BASIC FUNCTIONALITY}
-	%\begin{itemize}
-	%	\item xxx
-	%	\item xxx
+	\underline{BASIC FUNCTIONALITY}
+	\begin{itemize}
+		\item When the Adaptive Time Step (ATS) Package was activated, and failed time steps were rerun successfully using a shorter time step, the cumulative budgets reported in the listing files were incorrect.  Cumulative budgets included contributions from the failed time steps.  The program was fixed so that cumulative budgets are not tallied until the time step is finalized.
+		\item For multi-model groundwater flow simulations that use the Buoyancy (BUY) Package, variable-density terms were not included, in all cases, along the interface between two GWF models.  The program was corrected so that flows between two GWF models include the variable-density terms when the BUY Package is active.
 	%	\item xxx
-	%\end{itemize}
+	\end{itemize}
 
 	%\underline{INTERNAL FLOW PACKAGES}
 	%\begin{itemize}

src/Model/GroundWaterFlow/gwf-lak.f90

@@ -4305,15 +4305,15 @@ end subroutine lak_ot_dv
   !<
   subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     ! -- module
-    use TdisModule, only: totim
+    use TdisModule, only: totim, delt
     ! -- dummy
     class(LakType) :: this !< LakType object
     integer(I4B), intent(in) :: kstp !< time step number
     integer(I4B), intent(in) :: kper !< period number
     integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
     integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
     !
-    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim)
+    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
     !
     ! -- Return
     return

src/Model/GroundWaterFlow/gwf-maw.f90

@@ -2792,15 +2792,15 @@ end subroutine maw_ot_dv
   !<
   subroutine maw_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     ! -- module
-    use TdisModule, only: totim
+    use TdisModule, only: totim, delt
     ! -- dummy
     class(MawType) :: this !< MawType object
     integer(I4B), intent(in) :: kstp !< time step number
     integer(I4B), intent(in) :: kper !< period number
     integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
     integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
     !
-    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim)
+    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
     !
     ! -- Return
     return

src/Model/GroundWaterFlow/gwf-mvr.f90

@@ -638,6 +638,7 @@ subroutine mvr_ot_bdsummary(this, ibudfl)
     end do
     !
     ! -- Write the budget
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       call this%budget%budget_ot(kstp, kper, this%iout)
     end if

src/Model/GroundWaterFlow/gwf-sfr.f90

@@ -2591,15 +2591,15 @@ end subroutine sfr_ot_dv
   !<
   subroutine sfr_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     ! -- module
-    use TdisModule, only: totim
+    use TdisModule, only: totim, delt
     ! -- dummy
     class(SfrType) :: this !< SfrType object
     integer(I4B), intent(in) :: kstp !< time step number
     integer(I4B), intent(in) :: kper !< period number
     integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
     integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
     !
-    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim)
+    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
     !
     ! -- return
     return

src/Model/GroundWaterFlow/gwf-uzf.f90

@@ -1670,15 +1670,15 @@ end subroutine uzf_ot_dv
   !<
   subroutine uzf_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     ! -- module
-    use TdisModule, only: totim
+    use TdisModule, only: totim, delt
     ! -- dummy
     class(UzfType) :: this !< UzfType object
     integer(I4B), intent(in) :: kstp !< time step number
     integer(I4B), intent(in) :: kper !< period number
     integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
     integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
     !
-    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim)
+    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
     !
     ! -- Return
     return

src/Model/GroundWaterFlow/gwf.f90

@@ -1032,7 +1032,7 @@ subroutine gwf_ot_dv(this, idvsave, idvprint, ipflag)
   end subroutine gwf_ot_dv
 
   subroutine gwf_ot_bdsummary(this, ibudfl, ipflag)
-    use TdisModule, only: kstp, kper, totim
+    use TdisModule, only: kstp, kper, totim, delt
     class(GwfModelType) :: this
     integer(I4B), intent(in) :: ibudfl
     integer(I4B), intent(inout) :: ipflag
@@ -1052,6 +1052,7 @@ subroutine gwf_ot_bdsummary(this, ibudfl, ipflag)
     end if
 
     ! -- model budget summary
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       ipflag = 1
       call this%budget%budget_ot(kstp, kper, this%iout)

src/Model/GroundWaterTransport/gwt-ist.f90

@@ -634,6 +634,7 @@ subroutine ist_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     call this%budget%addentry(this%budterm, delt, budtxt, isuppress_output)
     !
     ! -- Write budget to list file
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       call this%budget%budget_ot(kstp, kper, iout)
     end if

src/Model/ParticleTracking/prt.f90

@@ -689,7 +689,7 @@ end subroutine prt_ot_dv
   !> @brief Print budget summary
   subroutine prt_ot_bdsummary(this, ibudfl, ipflag)
     ! -- modules
-    use TdisModule, only: kstp, kper, totim
+    use TdisModule, only: kstp, kper, totim, delt
     ! -- dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: ibudfl
@@ -705,6 +705,7 @@ subroutine prt_ot_bdsummary(this, ibudfl, ipflag)
     end do
 
     ! -- model budget summary
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       ipflag = 1
       ! -- model budget summary

src/Model/SurfaceWaterFlow/swf.f90

@@ -855,7 +855,7 @@ subroutine swf_ot_dv(this, idvsave, idvprint, ipflag)
   end subroutine swf_ot_dv
 
   subroutine swf_ot_bdsummary(this, ibudfl, ipflag)
-    use TdisModule, only: kstp, kper, totim
+    use TdisModule, only: kstp, kper, totim, delt
     class(SwfModelType) :: this
     integer(I4B), intent(in) :: ibudfl
     integer(I4B), intent(inout) :: ipflag
@@ -875,6 +875,7 @@ subroutine swf_ot_bdsummary(this, ibudfl, ipflag)
     ! end if
 
     ! -- model budget summary
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       ipflag = 1
       call this%budget%budget_ot(kstp, kper, this%iout)

src/Model/TransportModel/tsp-apt.f90

@@ -1067,15 +1067,15 @@ end subroutine apt_ot_dv
   !<
   subroutine apt_ot_bdsummary(this, kstp, kper, iout, ibudfl)
     ! -- module
-    use TdisModule, only: totim
+    use TdisModule, only: totim, delt
     ! -- dummy
     class(TspAptType) :: this !< TspAptType object
     integer(I4B), intent(in) :: kstp !< time step number
     integer(I4B), intent(in) :: kper !< period number
     integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
     integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
     !
-    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim)
+    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
     !
     ! -- Return
     return

src/Model/TransportModel/tsp-mvt.f90

@@ -484,7 +484,6 @@ subroutine mvt_ot_bdsummary(this, ibudfl)
     ! -- locals
     integer(I4B) :: i, j, n
     real(DP), allocatable, dimension(:) :: ratin, ratout
-! ------------------------------------------------------------------------------
     !
     ! -- Allocate and initialize ratin/ratout
     allocate (ratin(this%maxpackages), ratout(this%maxpackages))
@@ -525,6 +524,7 @@ subroutine mvt_ot_bdsummary(this, ibudfl)
     end do
     !
     ! -- Write the budget
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       call this%budget%budget_ot(kstp, kper, this%iout)
     end if

src/Model/TransportModel/tsp.f90

@@ -509,7 +509,7 @@ end subroutine tsp_ot_dv
   !! Loop through attached packages and write budget summaries
   !<
   subroutine tsp_ot_bdsummary(this, ibudfl, ipflag)
-    use TdisModule, only: kstp, kper, totim
+    use TdisModule, only: kstp, kper, totim, delt
     class(TransportModelType) :: this
     integer(I4B), intent(in) :: ibudfl
     integer(I4B), intent(inout) :: ipflag
@@ -528,6 +528,7 @@ subroutine tsp_ot_bdsummary(this, ibudfl, ipflag)
     end if
     !
     ! -- Model budget summary
+    call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       ipflag = 1
       call this%budget%budget_ot(kstp, kper, this%iout)

src/Utilities/Budget.f90

@@ -64,6 +64,7 @@ module BudgetModule
     procedure :: add_single_entry
     procedure :: add_multi_entry
     generic :: addentry => add_single_entry, add_multi_entry
+    procedure :: finalize_step
     procedure :: writecsv
     ! -- private
     procedure :: allocate_scalars
@@ -84,7 +85,6 @@ subroutine budget_cr(this, name_model)
     ! -- dummy
     type(BudgetType), pointer :: this !< BudgetType object
     character(len=*), intent(in) :: name_model !< name of the model
-! ------------------------------------------------------------------------------
     !
     ! -- Create the object
     allocate (this)
@@ -363,12 +363,11 @@ subroutine reset(this)
     integer(I4B) :: i
     !
     this%msum = 1
-    !
     do i = 1, this%maxsize
       this%vbvl(3, i) = DZERO
       this%vbvl(4, i) = DZERO
     end do
-    !
+
     ! -- Return
     return
   end subroutine reset
@@ -425,11 +424,6 @@ subroutine add_single_entry(this, rin, rout, delt, text, &
       end if
     end if
     !
-    if (iscv == 0) then
-      this%vbvl(1, this%msum) = this%vbvl(1, this%msum) + rin * delt
-      this%vbvl(2, this%msum) = this%vbvl(2, this%msum) + rout * delt
-    end if
-    !
     this%vbvl(3, this%msum) = rin
     this%vbvl(4, this%msum) = rout
     this%vbnm(this%msum) = adjustr(text)
@@ -499,11 +493,6 @@ subroutine add_multi_entry(this, budterm, delt, budtxt, &
         end if
       end if
       !
-      if (iscv == 0) then
-        this%vbvl(1, this%msum) = this%vbvl(1, this%msum) + budterm(1, i) * delt
-        this%vbvl(2, this%msum) = this%vbvl(2, this%msum) + budterm(2, i) * delt
-      end if
-      !
       this%vbvl(3, this%msum) = budterm(1, i)
       this%vbvl(4, this%msum) = budterm(2, i)
       this%vbnm(this%msum) = adjustr(budtxt(i))
@@ -524,6 +513,29 @@ subroutine add_multi_entry(this, budterm, delt, budtxt, &
     return
   end subroutine add_multi_entry
 
+  !> @ brief Update accumulators
+  !!
+  !! This must be called before any output is written
+  !! in order to update the accumulators in vbvl(1,:)
+  !! and vbl(2,:).
+  !<
+  subroutine finalize_step(this, delt)
+    ! -- modules
+    ! -- dummy
+    class(BudgetType) :: this !< BudgetType object
+    real(DP), intent(in) :: delt
+    ! -- local
+    integer(I4B) :: i
+    !
+    do i = 1, this%msum - 1
+      this%vbvl(1, i) = this%vbvl(1, i) + this%vbvl(3, i) * delt
+      this%vbvl(2, i) = this%vbvl(2, i) + this%vbvl(4, i) * delt
+    end do
+    !
+    ! -- Return
+    return
+  end subroutine finalize_step
+
   !> @ brief allocate scalar variables
   !!
   !!  Allocate scalar variables of this budget object