diff --git a/autotest/test_gwf_ats_lak01.py b/autotest/test_gwf_ats_lak01.py index 90d9a0ad55d..70dbac9a2fd 100644 --- a/autotest/test_gwf_ats_lak01.py +++ b/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)) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index fbab15b296b..4fd7e1b877d 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/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} diff --git a/src/Model/GroundWaterFlow/gwf-lak.f90 b/src/Model/GroundWaterFlow/gwf-lak.f90 index 4218ca8d3f2..e6386e6fcd5 100644 --- a/src/Model/GroundWaterFlow/gwf-lak.f90 +++ b/src/Model/GroundWaterFlow/gwf-lak.f90 @@ -4305,7 +4305,7 @@ 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 @@ -4313,7 +4313,7 @@ subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl) 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 diff --git a/src/Model/GroundWaterFlow/gwf-maw.f90 b/src/Model/GroundWaterFlow/gwf-maw.f90 index df40289b0d3..cf603f6ce42 100644 --- a/src/Model/GroundWaterFlow/gwf-maw.f90 +++ b/src/Model/GroundWaterFlow/gwf-maw.f90 @@ -2792,7 +2792,7 @@ 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 @@ -2800,7 +2800,7 @@ subroutine maw_ot_bdsummary(this, kstp, kper, iout, ibudfl) 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 diff --git a/src/Model/GroundWaterFlow/gwf-mvr.f90 b/src/Model/GroundWaterFlow/gwf-mvr.f90 index 0e9fd9594ce..8f35a23414f 100644 --- a/src/Model/GroundWaterFlow/gwf-mvr.f90 +++ b/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 diff --git a/src/Model/GroundWaterFlow/gwf-sfr.f90 b/src/Model/GroundWaterFlow/gwf-sfr.f90 index 27fc3e60f67..7c42b81742e 100644 --- a/src/Model/GroundWaterFlow/gwf-sfr.f90 +++ b/src/Model/GroundWaterFlow/gwf-sfr.f90 @@ -2591,7 +2591,7 @@ 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 @@ -2599,7 +2599,7 @@ subroutine sfr_ot_bdsummary(this, kstp, kper, iout, ibudfl) 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 diff --git a/src/Model/GroundWaterFlow/gwf-uzf.f90 b/src/Model/GroundWaterFlow/gwf-uzf.f90 index d8b7019fc6a..1ed322dcd92 100644 --- a/src/Model/GroundWaterFlow/gwf-uzf.f90 +++ b/src/Model/GroundWaterFlow/gwf-uzf.f90 @@ -1670,7 +1670,7 @@ 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 @@ -1678,7 +1678,7 @@ subroutine uzf_ot_bdsummary(this, kstp, kper, iout, ibudfl) 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 diff --git a/src/Model/GroundWaterFlow/gwf.f90 b/src/Model/GroundWaterFlow/gwf.f90 index 1caed0d0fec..5b7d6ee90c0 100644 --- a/src/Model/GroundWaterFlow/gwf.f90 +++ b/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) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index ad0a91ee4a0..88c975823b7 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/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 diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index ac3af4f9f4c..9dc4c1f1aee 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/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 diff --git a/src/Model/SurfaceWaterFlow/swf.f90 b/src/Model/SurfaceWaterFlow/swf.f90 index 75ade14ed92..d592696906d 100644 --- a/src/Model/SurfaceWaterFlow/swf.f90 +++ b/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) diff --git a/src/Model/TransportModel/tsp-apt.f90 b/src/Model/TransportModel/tsp-apt.f90 index b63a3b70569..47fa4611c35 100644 --- a/src/Model/TransportModel/tsp-apt.f90 +++ b/src/Model/TransportModel/tsp-apt.f90 @@ -1067,7 +1067,7 @@ 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 @@ -1075,7 +1075,7 @@ subroutine apt_ot_bdsummary(this, kstp, kper, iout, ibudfl) 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 diff --git a/src/Model/TransportModel/tsp-mvt.f90 b/src/Model/TransportModel/tsp-mvt.f90 index f16961666b2..8a52a1b50c2 100644 --- a/src/Model/TransportModel/tsp-mvt.f90 +++ b/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 diff --git a/src/Model/TransportModel/tsp.f90 b/src/Model/TransportModel/tsp.f90 index 15f997bff07..a3eb1f70589 100644 --- a/src/Model/TransportModel/tsp.f90 +++ b/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) diff --git a/src/Utilities/Budget.f90 b/src/Utilities/Budget.f90 index e3947909e80..b91f52d398a 100644 --- a/src/Utilities/Budget.f90 +++ b/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 diff --git a/src/Utilities/BudgetObject.f90 b/src/Utilities/BudgetObject.f90 index 87819b7dd10..83cced41090 100644 --- a/src/Utilities/BudgetObject.f90 +++ b/src/Utilities/BudgetObject.f90 @@ -459,7 +459,7 @@ end subroutine write_flowtable !> @brief Write the budget table !< - subroutine write_budtable(this, kstp, kper, iout, ibudfl, totim) + subroutine write_budtable(this, kstp, kper, iout, ibudfl, totim, delt) ! -- dummy class(BudgetObjectType) :: this integer(I4B), intent(in) :: kstp @@ -467,8 +467,10 @@ subroutine write_budtable(this, kstp, kper, iout, ibudfl, totim) integer(I4B), intent(in) :: iout integer(I4B), intent(in) :: ibudfl real(DP), intent(in) :: totim + real(DP), intent(in) :: delt ! ! -- Write the table + call this%budtable%finalize_step(delt) if (ibudfl /= 0) then call this%budtable%budget_ot(kstp, kper, iout) end if diff --git a/utils/zonebudget/src/zoneoutput.f90 b/utils/zonebudget/src/zoneoutput.f90 index 8e3c556d48b..6e304940632 100644 --- a/utils/zonebudget/src/zoneoutput.f90 +++ b/utils/zonebudget/src/zoneoutput.f90 @@ -216,6 +216,7 @@ subroutine zoneoutput_write(itime, kstp, kper, delt, totim, nbudterms, & call budobj(izone)%addentry(rin, rout, delt, txt) end if end do + call budobj(izone)%finalize_step(delt) call budobj(izone)%budget_ot(kstp, kper, iout) ! ! -- write line ending after each zone