From 9705a3077c92de94bf456447be4f3765b9e459ab Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Wed, 1 Nov 2023 07:09:12 -0500 Subject: [PATCH] ci: add simple gwfgwf exchange model with sfr mover (#1418) * ci: add simple gwfgwf exchange model with sfr mover * fix sfr gwtgwt test to add mover to exchange files * * modify simulation.py to add run_simulation keyword which allows a comparison simulation to be run but not compared --- autotest/simulation.py | 7 +- autotest/test_gwf_sfr01gwfgwf.py | 322 +++++++++++++++++++++++++++++++ autotest/test_gwt_sft01gwtgwt.py | 7 +- 3 files changed, 330 insertions(+), 6 deletions(-) create mode 100644 autotest/test_gwf_sfr01gwfgwf.py diff --git a/autotest/simulation.py b/autotest/simulation.py index 8471acd64e2..a4d8982ccdb 100644 --- a/autotest/simulation.py +++ b/autotest/simulation.py @@ -48,6 +48,7 @@ def __init__( api_func=None, mf6_regression=False, make_comparison=True, + run_comparison=None, simpath=None, ): msg = sfmt.format("Initializing test", name) @@ -65,6 +66,10 @@ def __init__( self.api_func = api_func self.mf6_regression = mf6_regression self.make_comparison = make_comparison + if run_comparison is None: + self.run_comparison = make_comparison + else: + self.run_comparison = run_comparison self.action = None # set htol for comparisons @@ -165,7 +170,7 @@ def setup(self, src, dst): def setup_comparison(self, src, dst, testModel=True): # evaluate if comparison should be made - if not self.make_comparison: + if not self.run_comparison: return # adjust htol if it is smaller than IMS outer_dvclose diff --git a/autotest/test_gwf_sfr01gwfgwf.py b/autotest/test_gwf_sfr01gwfgwf.py new file mode 100644 index 00000000000..5673751c4a7 --- /dev/null +++ b/autotest/test_gwf_sfr01gwfgwf.py @@ -0,0 +1,322 @@ +# Based on sft01 gwf model, but split into two gwf models test gwf-gwf and +# mvr. The single base model is split using the model splitter into two models. +# The single model is run as the regression model + +# The final split model look like: +# +# flow1 flow2 +# sfr 1 2 3 4 5 6 7 gwfgwf-mvr => 1 2 3 4 5 6 7 +# ------------- ------------- +# gwf 1 2 3 4 5 6 7 gwfgwf => 1 2 3 4 5 6 7 + + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + +# from flopy.mf6.utils import Mf6Splitter + +ex = ["sfr01gwfgwf"] + +# properties for single model combination +lx = 14.0 +lz = 1.0 +nlay = 1 +nrow = 1 +ncol = 14 +nper = 1 +delc = 1.0 +delr = lx / ncol +delz = lz / nlay +top = 0.0 +botm = [top - (k + 1) * delz for k in range(nlay)] +Kh = 20.0 +Kv = 20.0 + + +def build_simulation(idx, sim_ws, sim_type="single"): + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, + sim_ws=sim_ws, + ) + + tdis = flopy.mf6.ModflowTdis( + sim, + time_units="DAYS", + nper=nper, + ) + + # Flow solver + ims = flopy.mf6.ModflowIms( + sim, + complexity="simple", + print_option="ALL", + outer_dvclose=1e-6, + inner_dvclose=1e-6, + ) + + if sim_type == "single": + gwf_types = ("single",) + else: + gwf_types = ("left", "right") + for gwf_type in gwf_types: + gwf = build_gwf(sim, gwf_type=gwf_type) + + if sim_type != "single": + build_exchanges(sim) + + return sim + + +def build_gwf(sim, gwf_type="single"): + if gwf_type == "single": + nc = ncol + else: + nc = int(ncol / 2) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwf_type, + save_flows=True, + ) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=nc, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=0, + k=Kh, + k33=Kv, + ) + + # add chd to right edge + if gwf_type in ("single", "right"): + chdlist = [ + [(0, 0, nc - 1), 0.0], + ] + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdlist, + pname="chd_right", + ) + + # inject water into left edge + if gwf_type in ("single", "left"): + wellist = [ + [(0, 0, 0), 1.0], + ] + wel = flopy.mf6.ModflowGwfwel( + gwf, + stress_period_data=wellist, + pname="well_left", + ) + + # pak_data = [ [] []] + rlen = delr + rwid = delc + rgrd = 1.0 + rtp = 0.0 + rbth = 0.1 + rhk = 0.01 + rman = 1.0 + ustrf = 1.0 + ndv = 0 + pak_data = [] + for irno in range(nc): + ncon = 2 + if irno in [0, nc - 1]: + ncon = 1 + cellid = (0, 0, irno) + t = ( + irno, + cellid, + rlen, + rwid, + rgrd, + rtp, + rbth, + rhk, + rman, + ncon, + ustrf, + ndv, + ) + pak_data.append(t) + + con_data = [] + for irno in range(nc): + if irno == 0: + t = (irno, -(irno + 1)) + elif irno == nc - 1: + t = (irno, irno - 1) + else: + t = (irno, irno - 1, -(irno + 1)) + con_data.append(t) + + if gwf_type in ("single", "left"): + p_data = [ + (0, "INFLOW", 1.0), + ] + else: + p_data = None + + if gwf_type != "single": + mover = True + else: + mover = None + + sfr = flopy.mf6.modflow.ModflowGwfsfr( + gwf, + save_flows=True, + print_input=True, + print_flows=True, + print_stage=True, + mover=mover, + stage_filerecord=f"{gwf_type}.sfr.stg", + budget_filerecord=f"{gwf_type}.sfr.bud", + nreaches=nc, + packagedata=pak_data, + connectiondata=con_data, + perioddata=p_data, + pname=f"sfr_{gwf_type}", + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{gwf_type}.cbc", + head_filerecord=f"{gwf_type}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("BUDGET", "LAST")], + ) + + return gwf + + +def build_exchanges(sim): + # add a gwf-gwf exchange + gwfgwf_data = [ + ( + (0, 0, int(ncol / 2) - 1), + (0, 0, 0), + 1, + delr / 2.0, + delr / 2.0, + delc, + 0.0, + delr, + ) + ] + + # GWF-GWF + mvr_filerecord = "left-right.exg.mvr" + gwfgwf = flopy.mf6.ModflowGwfgwf( + sim, + exgtype="GWF6-GWF6", + nexg=len(gwfgwf_data), + exgmnamea="left", + exgmnameb="right", + exchangedata=gwfgwf_data, + auxiliary=["ANGLDEGX", "CDIST"], + dev_interfacemodel_on=False, + filename="left-right.exg", + ) + + # simulation GWF-GWF Mover + maxmvr, maxpackages = 1, 2 + mvrpack_sim = [["left", "sfr_left"], ["right", "sfr_right"]] + mvrspd = [ + [ + "left", + "sfr_left", + int(ncol / 2) - 1, + "right", + "sfr_right", + 0, + "FACTOR", + 1.00, + ] + ] + + gwfgwf.mvr.initialize( + modelnames=True, + maxmvr=maxmvr, + print_flows=True, + maxpackages=maxpackages, + packages=mvrpack_sim, + perioddata=mvrspd, + filename=mvr_filerecord, + ) + + +def build_model(idx, ws): + sim_ws = ws / "mf6" + sim_base = build_simulation(idx, sim_ws) + + sim = build_simulation(idx, ws, sim_type="split") + + return sim, sim_base + + +def eval_results(sim): + print("evaluating sfr stage results...") + + # base simulations stage + ws = sim.simpath + fpth = ws / "mf6/single.sfr.stg" + single_stage_obj = flopy.utils.HeadFile(fpth, text="STAGE") + single_stage = single_stage_obj.get_data().squeeze() + + stage = single_stage.copy() + + i1 = int(ncol / 2) + fpth = ws / "left.sfr.stg" + stage_obj = flopy.utils.HeadFile(fpth, text="STAGE") + v = stage_obj.get_data().squeeze() + stage[:i1] = v[:] + + fpth = ws / "right.sfr.stg" + stage_obj = flopy.utils.HeadFile(fpth, text="STAGE") + v = stage_obj.get_data().squeeze() + stage[i1:] = v[:] + + assert np.allclose(single_stage, stage), "sfr stages are not equal" + + +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = function_tmpdir + test = TestFramework() + test.build(build_model, idx, ws) + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + run_comparison=True, + ), + ws, + ) diff --git a/autotest/test_gwt_sft01gwtgwt.py b/autotest/test_gwt_sft01gwtgwt.py index ba63d36a3a2..f1651670e4e 100644 --- a/autotest/test_gwt_sft01gwtgwt.py +++ b/autotest/test_gwt_sft01gwtgwt.py @@ -60,7 +60,6 @@ def build_model(idx, ws): - name = "mf6sim" sim = flopy.mf6.MFSimulation( sim_name=name, @@ -180,8 +179,7 @@ def build_model(idx, ws): mvrspd = [ ["flow1", "sfr-1", ncol - 1, "flow2", "sfr-1", 0, "FACTOR", 1.00] ] - mvr = flopy.mf6.ModflowMvr( - sim, + gwfgwf.mvr.initialize( modelnames=True, maxmvr=maxmvr, print_flows=True, @@ -212,14 +210,13 @@ def build_model(idx, ws): # simulation GWT-GWT Mover if across_model_mvt_on: - mvt = flopy.mf6.modflow.ModflowGwtmvt(sim, filename=mvt_filerecord) + gwtgwt.mvt.initialize(filename=mvt_filerecord) regression = None return sim, regression def build_gwfgwt_combo(sim, gwfname, gwtname, icombo): - # create gwf model gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname)