From 35fe429625a5793a81bb7ae3614347ce90eb1aae Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D" Date: Tue, 2 Jan 2024 10:24:01 -0600 Subject: [PATCH 1/2] fix(gwfgwf): budget error when chd lives on border of exchange * close #1529 --- autotest/test_gwfgwf.py | 217 ++++++++++++++++++++++++++++++++ doc/ReleaseNotes/develop.tex | 8 +- src/Exchange/GwfGwfExchange.f90 | 87 ++++++++++++- 3 files changed, 302 insertions(+), 10 deletions(-) create mode 100644 autotest/test_gwfgwf.py diff --git a/autotest/test_gwfgwf.py b/autotest/test_gwfgwf.py new file mode 100644 index 00000000000..d2b36f6ea7c --- /dev/null +++ b/autotest/test_gwfgwf.py @@ -0,0 +1,217 @@ +""" +Test a constant head boundary assigned to one model as +also functioning as a constant head boundary assigned to +a connected model. The constant head term for model 1 +should be calculated to include the flow from model 2. +Also, the flowja budget term for the constant head cell +should have a correct residual of zero in the diagonal +position. + + 1 1 -1 gwf1 + - - - + 1 1 1 gwf2 + +We assert equality on qresidual being less than tolerance and +also that total constant head flow is correct. +""" + +import os +import pathlib as pl +import math +import flopy +import numpy as np +import pytest + +from framework import TestFramework + +cases = ["gwfgwf01", "gwfgwf01ifmod"] +ifmod = [False, True] + + +def build_models(idx, test): + sim = get_sim(idx, test.workspace) + return sim, None + + +def get_sim(idx, dir): + name = cases[idx] + + # solver data + nouter, ninner = 100, 300 + hclose, rclose, relax = 1.e-8, 1e-8, 0.97 + + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=dir + ) + + tdis = flopy.mf6.ModflowTdis( + sim, + time_units="DAYS", + nper=1, + perioddata=[(1.0, 1, 1)], + ) + + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + relaxation_factor=relax, + ) + + gwf0 = add_model(sim, "gwf1", top=0., add_chd=True, add_rch=True) + gwf1 = add_model(sim, "gwf2", top=-1., add_chd=False, add_rch=False) + gwfgwf = add_gwfexchange(sim, idx) + + return sim + + +def add_model(sim, modelname, top, add_chd, add_rch): + + nlay, nrow, ncol = 1, 1, 3 + botm = top - 1.0 + gwf = flopy.mf6.ModflowGwf(sim, modelname=modelname, save_flows=True) + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=1.0, + delc=1.0, + top=top, + botm=botm, + ) + ic = flopy.mf6.ModflowGwfic(gwf, strt=0.) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + save_flows=True, + icelltype=0, + k=1.0, + ) + + if add_chd: + chdlist = [(0, 0, ncol-1, 0.)] + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdlist) + + if add_rch: + rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.001) + + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{modelname}.hds", + budget_filerecord=f"{modelname}.cbc", + budgetcsv_filerecord=f"{modelname}.cbc.csv", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + ) + + return gwf + + +def add_gwfexchange(sim, idx): + ncol = 3 + delr = 1. + delc = 1. + dz = 1. + angldegx = 0.0 + cdist = 1.0 + gwfgwf_data = [ + [ + (0, 0, icol), + (0, 0, icol), + 0, + dz / 2.0, + dz / 2.0, + delr * delc, + angldegx, + cdist, + ] + for icol in range(ncol) + ] + gwfgwf = flopy.mf6.ModflowGwfgwf( + sim, + exgtype="GWF6-GWF6", + save_flows=True, + print_flows=True, + nexg=len(gwfgwf_data), + exgmnamea="gwf1", + exgmnameb="gwf2", + exchangedata=gwfgwf_data, + auxiliary=["ANGLDEGX", "CDIST"], + dev_interfacemodel_on=ifmod[idx], + ) + return gwfgwf + + +def check_output(idx, test): + sim = test.sims[0] + check_model(sim, 0) + check_model(sim, 1) + + +def check_model(sim, model_number): + print (f"Checking model output (gwf{model_number})") + gwf = sim.gwf[model_number] + sim_ws = pl.Path(sim.sim_path) + fpth = sim_ws / f"gwf{model_number + 1}.dis.grb" + grb = flopy.mf6.utils.MfGrdFile(fpth) + ia = grb.ia + + bobj = gwf.output.budget() + # print(bobj.list_records()) + budget_records = bobj.get_data(kstpkper=(0, 0)) + nrecords = len(budget_records) + print(f"detected {nrecords} items in budget file.") + for idx in range(nrecords): + budget_record = bobj.get_data(idx=idx)[0] + print(budget_record) + + model_budget = gwf.output.budgetcsv().data + pd = model_budget["PERCENT_DIFFERENCE"][0] + print("percent difference: ", model_budget["PERCENT_DIFFERENCE"]) + errmsg = f"Model percent difference is too large (pd)" + assert pd < 1.e-6, errmsg + + # check residual budget term in flowja diagonal position + fja = bobj.get_data(text="FLOW-JA-FACE")[0].flatten() + success = True + atol = 1.e-7 + for ipos in ia[:-1]: + print(ipos, fja[ipos]) + if fja[ipos] > atol: + success = False + assert success, f"flowja residual larger than tolerance ({atol})" + + # ensure that the constant head outflow is equal to the + # specified recharge + if model_number == 0: + chdflows = bobj.get_data(text="chd")[0] + vchd = -chdflows["q"][0] + vexgchd = model_budget["FLOW-JA-FACE-CHD(GWF-GWF_1)_OUT"][0] + v = vchd + vexgchd + errmsg = ( + f"Constant head outflow ({v}) is not equal to the " + f"specified recharge inflow (0.002)." + ) + assert math.isclose(v, 0.002), errmsg + + return + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 8b9bf02b093..fd7e5778419 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -55,9 +55,9 @@ % \item xxx %\end{itemize} - %\underline{EXCHANGES} - %\begin{itemize} - % \item xxx + \underline{EXCHANGES} + \begin{itemize} + \item A model budget error would occur when a constant-head cell in one model had a direct connection to an active cell in another model. For the model budget to be calculated correctly a new term called ``FLOW-JA-FACE-CHD'' was added to the GWF model budget. This term is only included in the table when the GWF Model is connected to another GWF Model using a GWF-GWF Exchange. % \item xxx % \item xxx - %\end{itemize} + \end{itemize} diff --git a/src/Exchange/GwfGwfExchange.f90 b/src/Exchange/GwfGwfExchange.f90 index 380e46d3ce3..2cc6e9c0281 100644 --- a/src/Exchange/GwfGwfExchange.f90 +++ b/src/Exchange/GwfGwfExchange.f90 @@ -99,6 +99,7 @@ module GwfGwfExchangeModule procedure, private :: condcalc procedure, private :: rewet procedure, private :: qcalc + procedure, private :: gwf_gwf_chd_bd procedure :: gwf_gwf_bdsav procedure, private :: gwf_gwf_bdsav_model procedure, private :: gwf_gwf_df_obs @@ -755,17 +756,21 @@ subroutine gwf_gwf_add_to_flowja(this) do i = 1, this%nexg ! if (associated(this%gwfmodel1)) then - flow = this%simvals(i) n = this%nodem1(i) - idiag = this%gwfmodel1%ia(n) - this%gwfmodel1%flowja(idiag) = this%gwfmodel1%flowja(idiag) + flow + if (this%gwfmodel1%ibound(n) > 0) then + flow = this%simvals(i) + idiag = this%gwfmodel1%ia(n) + this%gwfmodel1%flowja(idiag) = this%gwfmodel1%flowja(idiag) + flow + end if end if ! if (associated(this%gwfmodel2)) then - flow = -this%simvals(i) n = this%nodem2(i) - idiag = this%gwfmodel2%ia(n) - this%gwfmodel2%flowja(idiag) = this%gwfmodel2%flowja(idiag) + flow + if (this%gwfmodel2%ibound(n) > 0) then + flow = -this%simvals(i) + idiag = this%gwfmodel2%ia(n) + this%gwfmodel2%flowja(idiag) = this%gwfmodel2%flowja(idiag) + flow + end if end if ! end do @@ -928,6 +933,10 @@ subroutine gwf_gwf_bd(this, icnvg, isuppress_output, isolnid) call this%gwfmodel2%model_bdentry(budterm, budtxt, this%name) end if ! + ! -- Add any flows from one model into a constant head in another model + ! as a separate budget term called FLOW-JA-FACE-CHD + call this%gwf_gwf_chd_bd() + ! ! -- Call mvr bd routine if (this%inmvr > 0) call this%mvr%mvr_bd() ! @@ -935,6 +944,72 @@ subroutine gwf_gwf_bd(this, icnvg, isuppress_output, isolnid) return end subroutine gwf_gwf_bd + !> @ brief gwf-gwf-chd-bd + !! + !! Account for flow from an external model into a chd cell + !< + subroutine gwf_gwf_chd_bd(this) + ! -- modules + use ConstantsModule, only: DZERO, LENBUDTXT, LENPACKAGENAME + use BudgetModule, only: rate_accumulator + ! -- dummy + class(GwfExchangeType) :: this !< GwfExchangeType + ! -- local + character(len=LENBUDTXT), dimension(1) :: budtxt + integer(I4B) :: n + integer(I4B) :: i + real(DP), dimension(2, 1) :: budterm + real(DP) :: ratin, ratout + real(DP) :: q + ! + ! -- initialize + budtxt(1) = 'FLOW-JA-FACE-CHD' + ! + ! -- Add the constant-head budget terms for flow from model 2 into model 1 + if (associated(this%gwfmodel1)) then + ratin = DZERO + ratout = DZERO + do i = 1, this%nexg + n = this%nodem1(i) + if (this%gwfmodel1%ibound(n) < 0) then + q = this%simvals(i) + if (q > DZERO) then + ratout = ratout + q + else + ratin = ratin - q + end if + end if + end do + budterm(1, 1) = ratin + budterm(2, 1) = ratout + call this%gwfmodel1%model_bdentry(budterm, budtxt, this%name) + end if + ! + ! -- Add the constant-head budget terms for flow from model 1 into model 2 + if (associated(this%gwfmodel2)) then + ratin = DZERO + ratout = DZERO + do i = 1, this%nexg + n = this%nodem2(i) + if (this%gwfmodel2%ibound(n) < 0) then + ! -- flip flow sign as flow is relative to model 1 + q = - this%simvals(i) + if (q > DZERO) then + ratout = ratout + q + else + ratin = ratin - q + end if + end if + end do + budterm(1, 1) = ratin + budterm(2, 1) = ratout + call this%gwfmodel2%model_bdentry(budterm, budtxt, this%name) + end if + ! + ! -- Return + return + end subroutine gwf_gwf_chd_bd + !> @ brief Budget save !! !! Output individual flows to listing file and binary budget files From f8269e3238e1f7606712bc4e58342d45f067abad Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D" Date: Tue, 2 Jan 2024 10:29:55 -0600 Subject: [PATCH 2/2] fprettify --- src/Exchange/GwfGwfExchange.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Exchange/GwfGwfExchange.f90 b/src/Exchange/GwfGwfExchange.f90 index 2cc6e9c0281..e48b7be02dc 100644 --- a/src/Exchange/GwfGwfExchange.f90 +++ b/src/Exchange/GwfGwfExchange.f90 @@ -993,7 +993,7 @@ subroutine gwf_gwf_chd_bd(this) n = this%nodem2(i) if (this%gwfmodel2%ibound(n) < 0) then ! -- flip flow sign as flow is relative to model 1 - q = - this%simvals(i) + q = -this%simvals(i) if (q > DZERO) then ratout = ratout + q else