From c80352996f90aa29801a8e21fd736e59d572934a Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Thu, 16 May 2024 12:42:47 -0700 Subject: [PATCH 1/8] fix(gwt-sft.f90): ensure that the size of an array is > 0 before attempting to deallocate --- autotest/test_gwe_drycell_cnd3.py | 516 ++++++++++++++++++++++ autotest/test_gwf_sfrsft_gwdischrg.py | 266 +++++++++++ autotest/test_gwt_uztmvt2x1.py | 607 ++++++++++++++++++++++++++ src/Utilities/BudgetTerm.f90 | 12 +- 4 files changed, 1396 insertions(+), 5 deletions(-) create mode 100644 autotest/test_gwe_drycell_cnd3.py create mode 100644 autotest/test_gwf_sfrsft_gwdischrg.py create mode 100644 autotest/test_gwt_uztmvt2x1.py diff --git a/autotest/test_gwe_drycell_cnd3.py b/autotest/test_gwe_drycell_cnd3.py new file mode 100644 index 00000000000..0a72cedb324 --- /dev/null +++ b/autotest/test_gwe_drycell_cnd3.py @@ -0,0 +1,516 @@ +""" +Test problem for GWE: Specifically, this tests conduction between an APT feature +hosted in a drycell with and without NEWTON activated + + Model configuration: + + ~: Represents + conduction + +----//----+ + / ~//~ /| + / ~//~ / | + / ~//~ / | + +----//----+ + + Dry cell -> | // | /| + | | / | + | |/ | + +----------+ + + water table | | / + | | | / + -----+----------+--- + +----------+ + +""" + +# Imports + +import os +import numpy as np +import pytest +import flopy + +from framework import TestFramework + + +# Monotonicity function +def isMonotonic(A): + x, y = [], [] + x.extend(A) + y.extend(A) + x.sort() + y.sort(reverse=True) + if np.all(x == A) or np.all(y == A): + return True + return False + + +# Base simulation and model name and workspace + +scheme = "UPSTREAM" +# scheme = "TVD" + +cases = ["sfecnd"] + + +# Model units +length_units = "meters" +time_units = "days" + +# Table MODFLOW 6 GWE comparison + +nrow = 1 +ncol = 1 +nlay = 2 +top = 2 +bot = np.array([[[1.0]], [[0.0]]], dtype=float) +strthd = 0.1 # Starting head ($m$) +delr = 1.0 # Column width ($m$) +delc = 1.0 # Row width ($m$) +k11 = 1.0 # Horizontal hydraulic conductivity ($m/d$) +ss = 1e-6 # Specific storage +sy = 0.20 # Specific Yield +prsity = 0.20 # Porosity +nper = 4 # Number of periods +perlen = [1, 1000, 1, 1000] # Simulation time ($days$) +nstp = [1, 10, 1, 10] # 10 day transient time steps +steady = {0: False} +transient = {0: True} + +# sfr data +nreaches = 1 +rlen = 1.0 +rwid = 0.1 +roughness = 0.01 +rbth = 0.1 +rhk = 0.0 +slope = 0.001 +ustrf = 1.0 +ndv = 0 +nconn = 0 +cellid = (0, 0, 0) +sfr_packagedata = [] +rp = [ + 0, + cellid, + rlen, + rwid, + slope, + top, + rbth, + rhk, + roughness, + nconn, + ustrf, + ndv, +] +sfr_packagedata.append(rp) + +# sfe data +strm_temp = 18.0 # ($C$) +K_therm_strmbed = 2.0 # ($W/m/C$) +rbthcnd = 0.0001 # ($m$) + +# Set some static model parameter values + +k33 = k11 # Vertical hydraulic conductivity ($m/d$) +idomain = 1 # All cells included in the simulation +iconvert = 1 # All cells are convertible + +icelltype = 1 # Cell conversion type (>1: unconfined) + +# GWE related parameters +strt_temp = 4.0 +dispersivity = 0.0 # dispersion (remember, 1D model) +ktw = 0.5918 +kts = 0.2700 +rhos = 1500.0 +rhow = 1000.0 +cps = 760.0 +cpw = 4183.0 +lhv = 2454.0 + +# Set solver parameter values (and related) +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-10, 1.0 +ttsmult = 1.0 + +# Set up temporal data used by TDIS file +tdis_rc = [] +for i in np.arange(nper): + tdis_rc.append((perlen[i], nstp[i], ttsmult)) + +# ### Create MODFLOW 6 GWE MT3DMS Example 1 Boundary Conditions +# +# No GWF, only Heat conduction simulated + + +def add_gwf_model(sim, gwfname, newton=False): + + # Instantiating MODFLOW 6 groundwater flow model + if newton: + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + newtonoptions="NEWTON", + save_flows=True, + model_nam_file="{}.nam".format(gwfname), + ) + else: + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + save_flows=True, + model_nam_file="{}.nam".format(gwfname), + ) + + # Instantiating MODFLOW 6 solver for flow model + imsgwf = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwfname), + ) + sim.register_ims_package(imsgwf, [gwfname]) + + # Instantiating MODFLOW 6 discretization package + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=bot, + idomain=1, + pname="DIS", + filename="{}.dis".format(gwfname), + ) + + # Instantiating MODFLOW 6 storage package + flopy.mf6.ModflowGwfsto( + gwf, + ss=ss, + sy=sy, + iconvert=iconvert, + steady_state=steady, + transient=transient, + pname="STO", + filename="{}.sto".format(gwfname), + ) + + # Instantiating MODFLOW 6 node-property flow package + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=icelltype, + k=k11, + k33=k33, + save_specific_discharge=True, + pname="NPF", + filename="{}.npf".format(gwfname), + ) + + # Instantiating MODFLOW 6 initial conditions package for flow model + flopy.mf6.ModflowGwfic( + gwf, + strt=strthd, + filename="{}.ic".format(gwfname), + ) + + # Instantiating MODFLOW 6 output control package for flow model + flopy.mf6.ModflowGwfoc( + gwf, + pname="OC", + head_filerecord="{}.hds".format(gwfname), + budget_filerecord="{}.cbc".format(gwfname), + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # sfr file + connectiondata = [] + rc = [0] + connectiondata.append(rc) + + sfr_perioddata = [ + [0, "inflow", 1.0], + ] + + pname = "SFR-" + gwfname[-1] + sfr = flopy.mf6.ModflowGwfsfr( + gwf, + print_stage=True, + print_flows=True, + print_input=True, + save_flows=True, + nreaches=nreaches, + packagedata=sfr_packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + pname=pname, + filename="{}.sfr".format(gwfname), + ) + fname = f"{gwfname}.sfr.obs" + sfr_obs = { + f"{fname}.sfrobs": [ + ("inflow", "ext-inflow", 1), + ("outflow", "ext-outflow", 1), + ] + } + sfr.obs.initialize(filename=fname, print_input=True, continuous=sfr_obs) + + return sim + + +def add_gwe_model(sim, gwename): + + gwe = flopy.mf6.ModflowGwe( + sim, modelname=gwename, model_nam_file="{}.nam".format(gwename) + ) + gwe.name_file.save_flows = True + + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwename), + ) + sim.register_ims_package(imsgwe, [gwe.name]) + + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwedis( + gwe, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=bot, + idomain=1, + pname="DIS", + filename="{}.dis".format(gwename), + ) + + # Instantiating MODFLOW 6 transport initial concentrations + flopy.mf6.ModflowGweic( + gwe, strt=strt_temp, pname="IC", filename="{}.ic".format(gwename) + ) + + # Instantiating MODFLOW 6 transport advection package + flopy.mf6.ModflowGweadv( + gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename) + ) + + # Instantiating MODFLOW 6 transport dispersion package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + alh=dispersivity, + ath1=dispersivity, + ktw=ktw * 86400, + kts=kts * 86400, + pname="CND", + filename="{}.cnd".format(gwename), + ) + + # Instantiating MODFLOW 6 transport mass storage package (formerly "reaction" package in MT3DMS) + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + cps=cps, + rhos=rhos, + packagedata=[cpw, rhow, lhv], + pname="EST", + filename="{}.est".format(gwename), + ) + + # Instantiating MODFLOW 6 source/sink mixing package for dealing with + # auxiliary temperature specified in WEL boundary package. + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm( + gwe, + sources=sourcerecarray, + pname="SSM", + filename="{}.ssm".format(gwename), + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGweoc( + gwe, + pname="OC", + budget_filerecord="{}.cbc".format(gwename), + temperature_filerecord="{}.ucn".format(gwename), + temperatureprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + ) + + # Instantiate Streamflow Energy Transport package + sfe_packagedata = [] + t = (0, strm_temp, K_therm_strmbed, rbthcnd) + sfe_packagedata.append(t) + + sfe_perioddata = [] + sfe_perioddata.append((0, "INFLOW", strm_temp)) + flwpckname = "SFR-" + gwename[-1] + + flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=True, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfe_packagedata, + reachperioddata=sfe_perioddata, + flow_package_name=flwpckname, + pname="SFE", + filename="{}.sfe".format(gwename), + ) + + return sim + + +def build_models(idx, test): + + # Base MF6 GWF model type + ws = test.workspace + name = cases[idx] + + print("Building MF6 model...()".format(name)) + + # generate names for each model + gwfname1 = "gwf-" + name + "nwt1" + gwfname2 = "gwf-" + name + "non2" + gwename1 = "gwe-" + name + "nwt1" + gwename2 = "gwe-" + name + "non2" + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + # Instantiating MODFLOW 6 time discretization + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_rc, time_units=time_units + ) + + # Build two flow models, one with NWT, one without + sim = add_gwf_model(sim, gwfname1, newton=True) + sim = add_gwf_model(sim, gwfname2, newton=False) + + # Add GWE models for each of the flow models above + sim = add_gwe_model(sim, gwename1) + sim = add_gwe_model(sim, gwename2) + + # Add the flow-transport exchanges + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname1, + exgmnameb=gwename1, + pname="GWFGWE1", + filename="{}.gwfgwe1".format(gwename1), + ) + + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname2, + exgmnameb=gwename2, + pname="GWFGWE2", + filename="{}.gwfgwe2".format(gwename2), + ) + + return sim, None + + +def check_output(idx, test): + print("check results...") + + # confirm that a single reach model is working as expected + # test_gwe_drycell_cnd4 test ensures that multi-reach setup is working as well + # read transport results from GWE model + name = cases[idx] + gwename1 = "gwe-" + name + "nwt1" + gwename2 = "gwe-" + name + "non2" + + fpth_nwt = os.path.join(test.workspace, f"{gwename1}.ucn") + fpth_non = os.path.join(test.workspace, f"{gwename2}.ucn") + + try: + # load temperatures + tobj_nwt = flopy.utils.HeadFile( + fpth_nwt, precision="double", text="TEMPERATURE" + ) + temp_nwt = tobj_nwt.get_alldata() + except: + assert False, f'could not load temperature data from "{fpth_nwt}"' + + try: + # load temperatures + tobj_non = flopy.utils.HeadFile( + fpth_non, precision="double", text="TEMPERATURE" + ) + temp_non = tobj_non.get_alldata() + except: + assert False, f'could not load temperature data from "{fpth_non}"' + + # Ensure constant temperatures are initiated properly in the 1st and 3rd + # stress periods, which are separated by period of "turning off" the + # constant temperature boundary + msg0 = ( + "Grid cell temperatures for the Newton and non-Newton base are " + "different and should NOT be" + ) + assert np.allclose(temp_nwt, temp_non), msg0 + + msg1 = ( + "All layer 2 cells should be strictly colder than the cell above " + "owing to a conductive flux of energy from the stream into the " + "upper layer, but that is not the case" + ) + assert np.all(np.diff(temp_nwt.squeeze(), axis=1) < 0), msg1 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(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/autotest/test_gwf_sfrsft_gwdischrg.py b/autotest/test_gwf_sfrsft_gwdischrg.py new file mode 100644 index 00000000000..d905697b292 --- /dev/null +++ b/autotest/test_gwf_sfrsft_gwdischrg.py @@ -0,0 +1,266 @@ +# Test groundwater discharge to a stream and then go on to test +# that a transport model with a single reach works. + +import math +import pathlib as pl + +import flopy +import numpy as np +import pytest + +from framework import TestFramework + +cases = ["sfr-gwfout", "sfr-gwf-trnsprt"] + +length_units = "m" +time_units = "sec" + +nrow = 1 +ncol = 1 +nlay = 1 +delr = delc = 1.0 + +nouter, ninner = 100, 300 +dvclose, rclose, relax = 1e-10, 1e-10, 1.0 + + +def add_gwt_model(sim, gwtname): + gwt = flopy.mf6.ModflowGwt( + sim, modelname=gwtname, model_nam_file="{}.nam".format(gwtname) + ) + gwt.name_file.save_flows = True + + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=dvclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=dvclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwt.name]) + + flopy.mf6.ModflowGwtdis( + gwt, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0.0, + botm=-100.0, + ) + + flopy.mf6.ModflowGwtic(gwt, strt=1.0) + flopy.mf6.ModflowGwtmst(gwt, porosity=0.2, filename=f"{gwtname}.mst") + + sourcerecarray = [("GHB", "AUX", "CONCENTRATION")] + flopy.mf6.ModflowGwtssm( + gwt, + sources=sourcerecarray, + pname="SSM", + filename="{}.ssm".format(gwtname), + ) + + # Instantiate Streamflow Transport package + sft_packagedata = [] + t = (0, 1.0) + sft_packagedata.append(t) + + sft_perioddata = [] + sft_perioddata.append((0, "INFLOW", 0.0)) + flwpckname = "SFR" + + flopy.mf6.modflow.ModflowGwtsft( + gwt, + boundnames=False, + save_flows=True, + print_input=True, + print_flows=True, + print_concentration=True, + concentration_filerecord=gwtname + ".sft.bin", + budget_filerecord=gwtname + ".sft.bud", + packagedata=sft_packagedata, + reachperioddata=sft_perioddata, + flow_package_name=flwpckname, + pname="SFT", + filename="{}.sft".format(gwtname), + ) + + flopy.mf6.ModflowGwtoc( + gwt, + budget_filerecord=f"{gwtname}.cbc", + concentration_filerecord=f"{gwtname}.ucn", + concentrationprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + ) + + return sim + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + flopy.mf6.ModflowTdis(sim, time_units=time_units) + flopy.mf6.ModflowIms( + sim, + inner_dvclose=1e-5, + inner_hclose=1e-6, + ) + + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name, + ) + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0.0, + botm=-100.0, + ) + flopy.mf6.ModflowGwfnpf( + gwf, + icelltype=1, # >0 means saturated thickness varies with computed head + ) + flopy.mf6.ModflowGwfic(gwf, strt=1.0) + flopy.mf6.ModflowGwfghb( + gwf, + auxiliary="CONCENTRATION", + stress_period_data=[((0, 0, 0), 1.0, 1e6, 1.0)], + pname="GHB", + ) + + # sfr data + # + package_data = [ + (0, (0, 0, 0), delr, 1.0, 1e-3, 0.0, 1.0, 1.0, 0.001, 0, 0.0, 0) + ] + connection_data = [(0)] + + sfr_obs = { + f"{name}.sfr.obs.csv": [ + ("gwf", "sfr", (0,)), + ("outflow", "ext-outflow", (0,)), + ("depth", "depth", (0,)), + ], + "filename": name + ".sfr.obs", + } + + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + nreaches=1, + packagedata=package_data, + connectiondata=connection_data, + observations=sfr_obs, + pname="SFR", + ) + + if idx > 0: + gwtname = "gwt-sft" + sim = add_gwt_model(sim, gwtname) + + # Add the flow-transport exchanges + flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=name, + exgmnameb=gwtname, + pname="GWFGWT1", + filename="{}.gwfgwt1".format(gwtname), + ) + + return sim, None + + +def check_output(idx, test): + answer = np.array( + [ + 1.0, + -0.92094535738673577, + -0.92094535738673577, + 0.79053721667952215e-1, + ] + ) + obs_pth = pl.Path(f"{test.workspace}/{cases[idx]}.sfr.obs.csv") + sim_data = flopy.utils.Mf6Obs(obs_pth).get_data() + data_names = sim_data.dtype.names + for ct, name in enumerate(data_names): + assert np.allclose( + sim_data[name][0], answer[ct] + ), f"simulated sfr {name} results do not match answer" + + if idx > 0: + sft_ans = np.array([[[[1.0]]]]) + + gwtname = "gwt-sft" + sft_obs_fl = pl.Path(f"{test.workspace}/{gwtname}.sft.bin") + try: + # load simulated concentration in SFT + cobj = flopy.utils.HeadFile( + sft_obs_fl, text="CONCENTRATION" # precision="double" + ) + sim_conc_sft = cobj.get_alldata() + except: + assert ( + False + ), f'could not load concentration data from "{sft_obs_fl}"' + + gwt_sim_conc = pl.Path(f"{test.workspace}/{gwtname}.ucn") + try: + # load simulated concentration of groundwater + cobj = flopy.utils.HeadFile( + gwt_sim_conc, text="CONCENTRATION" # precision="double" + ) + conc_gw = cobj.get_alldata() + except: + assert ( + False + ), f'could not load temperature data from "{sft_obs_fl}"' + + msg0 = "The simulation is not matching the established answer" + msg1 = ( + "Groundwater discharge is the only source of flow in the " + "channel. Thus, the gw and sft water should have the same " + "concentration, but don't" + ) + assert np.allclose(sft_ans, sim_conc_sft), msg0 + assert np.allclose(conc_gw, sim_conc_sft), msg1 + + +@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/autotest/test_gwt_uztmvt2x1.py b/autotest/test_gwt_uztmvt2x1.py new file mode 100644 index 00000000000..0b78660ebfa --- /dev/null +++ b/autotest/test_gwt_uztmvt2x1.py @@ -0,0 +1,607 @@ +""" +2x1 test problem for GWE + +Test the rejected infiltration mover transport values in a sim with only a +single reach + +Model configuration: 1 SFR reaches exist in the first row. + 1 UZF objects exist in the second row + 1 MVR/MVT connections transfer rejected infiltration + from UZF to SFR + + + Row 1 +----------+ + SFR / /| +Channel-> ==================== + / ^ / | Profile View : + / / / | ------------ + +----+-----+ | uzf rej inf. --> mvt --> sfr/sft + / / /| + +-----------+ + Row 2 / rej. / | / | +---- ----+ + / inf. / | / | | \_/ | + / / | / | | | + +----------+ |/ | | | + | | + +-----------+ | + | | / +-----------+ + | | / + | | / + | |/ + +----------+ + +""" + +# Imports + +import os +import numpy as np +import pytest +import flopy + +from framework import TestFramework + +# Base simulation and model name and workspace + +scheme = "UPSTREAM" +# scheme = "TVD" + +cases = [ + "uztmvt" +] # 2-cell model, horizontally connected with staggered alignment + +nrow = 2 +ncol = 1 +nlay = 1 +top = np.array([[[1.0], [1.5]]], dtype=float) +bot = np.array([[[0.0], [0.5]]], dtype=float) +strthd = np.array([[[0.01], [0.51]]], dtype=float) + +# Model units +length_units = "meters" +time_units = "days" + +# Table MODFLOW 6 GWE comparison to MT3DMS + +delr = 1.0 # Column width ($m$) +delc = 1.0 # Row width ($m$) +k11 = 1.0e3 # Horizontal hydraulic conductivity ($m/d$) +vks = 0.01 # vertical hydraulic conductivity of the unsaturated zone +ss = 1e-6 # Specific storage +sy = 0.30 # Specific Yield +prsity = 0.20 # Porosity +nper = 2 # Number of periods +perlen = [1, 1] # Simulation time ($days$) +nstp = [1, 1] # 10 day transient time steps +ttsmult = 1.0 +steady = {0: False} +transient = {0: True} + +# Set some static model parameter values + +k33 = k11 # Vertical hydraulic conductivity ($m/d$) +idomain = 1 # All cells included in the simulation +iconvert = 1 # All cells are convertible + +icelltype = 1 # Cell conversion type (>1: unconfined) + +chdlist = [] +chdlist.append([(0, 0, 0), 0.51]) + +# Set some static transport related model parameter values +strt_conc = 10.0 +strt_uz_conc = 1.0 +dispersivity = 0.0 + +# UZF related parameters +thtr = 0.05 # Residual water content +thts = sy # Saturated water content +thti = 0.05 # Initial water content (unsaturated zone) +eps = 7.1 # Brooks-Corey epsilon + +# GWT related parameters +al = 0.0 # Longitudinal dispersivity ($m$) +rhob = 1.5 # Bulk density of layer 1 ($g/cm^3$) +Kd = 0.176 # Distribution coefficient ($cm^3/g$) + +# UZF/UZT related input +surfdep = 1.0e-5 +finf1 = 1.01 +finf2 = 2.01 +pet = 0.0 +extdp = 0.5 + +uzf_pkdat = [[0, (0, 1, 0), 1, -1, surfdep, vks, thtr, thts, thti, eps]] +uzf_spd = { + 0: [ + [0, finf1, pet, extdp, thtr, 0.0, 0.0, 0.0], + ], + 1: [[0, finf2, pet, extdp, thtr, 0.0, 0.0, 0.0]], +} +concCell = [] +concCell.append(11.1) + +uzt_pkdat = [(0, strt_uz_conc)] # ifno, strt_conc +uzt_perdat = [(0, "INFILTRATION", concCell[0])] + +# SFR/SFT related input +conn_dat = [0] +sfr_pkdat = [] +nreaches = 1 +rlen = 1.0 +rwid = 1.0 +roughness = 0.03 +rbth = 0.1 +rhk = 0.0 +slope = 0.001 +ustrf = 1.0 +ndv = 0 +rp = [ + 0, + (0, 0, 0), + rlen, + rwid, + slope, + top[0, 0, 0], + rbth, + rhk, + roughness, + 0, # ncon + ustrf, + 0, # ndv +] +sfr_pkdat.append(rp) +sfr_perdat = {0: [0, "INFLOW", 1.0]} +sft_pkdat = [0, 0.0] +sft_perdat = {0: [0, "INFLOW", 0.0]} + +# MVR input +mvr_pkdat = [] +mvr_pkdat.append(["UZF-1", 0, "SFR-1", 0, "FACTOR", 1.0]) +mvr_packages = [ + ("UZF-1",), + ("SFR-1",), + # ("DRN-1",), +] + +# Set solver parameter values (and related) +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-10, 1.0 + +# Set up temporal data used by TDIS file +tdis_rc = [] +for i in np.arange(nper): + tdis_rc.append((perlen[i], nstp[i], ttsmult)) + +# ### Create MODFLOW 6 GWT +# +# No GWF, only Heat conduction simulated + + +def build_models(idx, test): + + # Base MF6 GWF model type + ws = test.workspace + name = cases[idx] + + print("Building MF6 model...()".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + # Instantiating MODFLOW 6 time discretization + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_rc, time_units=time_units + ) + + # Instantiating MODFLOW 6 groundwater flow model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + save_flows=True, + model_nam_file="{}.nam".format(gwfname), + ) + + # Instantiating MODFLOW 6 solver for flow model + imsgwf = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="CG", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwfname), + ) + sim.register_ims_package(imsgwf, [gwfname]) + + # Instantiating MODFLOW 6 discretization package + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=bot, + idomain=1, + pname="DIS-1", + filename="{}.dis".format(gwfname), + ) + + # Instantiating MODFLOW 6 storage package + # flopy.mf6.ModflowGwfsto( + # gwf, + # ss=ss, + # sy=sy, + # iconvert=iconvert, + # steady_state=steady, + # transient=transient, + # pname="STO", + # filename="{}.sto".format(gwfname), + # ) + + # CHD + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdlist, + pname="CHD", + filename="{}.chd".format(gwfname), + ) + + # Instantiating MODFLOW 6 node-property flow package + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=icelltype, + k=k11, + k33=k33, + save_specific_discharge=True, + pname="NPF-1", + filename="{}.npf".format(gwfname), + ) + + flopy.mf6.ModflowGwfsto( + gwf, + ss=ss, + sy=sy, + transient={0: True}, + pname="STO", + filename="{}.sto".format(gwfname), + ) + + # Instantiating MODFLOW 6 initial conditions package for flow model + flopy.mf6.ModflowGwfic( + gwf, + strt=strthd, + filename="{}.ic".format(gwfname), + ) + + # UZF + flopy.mf6.ModflowGwfuzf( + gwf, + mover=True, + print_input=True, + print_flows=True, + save_flows=True, + boundnames=False, + simulate_et=False, + unsat_etwc=False, + ntrailwaves=15, + nwavesets=40, + nuzfcells=len(uzf_pkdat), + packagedata=uzf_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{gwfname}.uzf.bud", + pname="UZF-1", + filename=f"{gwfname}.uzf", + ) + + # SFR + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=86400.0, + budget_filerecord=f"{gwfname}.sfr.bud", + mover=True, + nreaches=nreaches, + packagedata=sfr_pkdat, + connectiondata=conn_dat, + perioddata=sfr_perdat, + pname="SFR-1", + filename="{}.sfr".format(gwfname), + ) + + # MVR + flopy.mf6.ModflowGwfmvr( + gwf, + maxmvr=len(mvr_pkdat), + budget_filerecord=f"{gwfname}.mvr.bud", + maxpackages=len(mvr_packages), + print_flows=True, + packages=mvr_packages, + perioddata=mvr_pkdat, + ) + + # Instantiating MODFLOW 6 output control package for flow model + flopy.mf6.ModflowGwfoc( + gwf, + pname="OC-1", + head_filerecord="{}.hds".format(gwfname), + budget_filerecord="{}.cbc".format(gwfname), + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # ---------------------------------------------------- + # Instantiating MODFLOW 6 GWT model + # ---------------------------------------------------- + gwt = flopy.mf6.ModflowGwt( + sim, modelname=gwtname, model_nam_file="{}.nam".format(gwtname) + ) + gwt.name_file.save_flows = True + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwt.name]) + + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwtdis( + gwt, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=bot, + idomain=1, + pname="DIS-1", + filename="{}.dis".format(gwtname), + ) + + # Instantiating MODFLOW 6 transport initial concentrations + flopy.mf6.ModflowGwtic( + gwt, strt=strt_conc, pname="IC-1", filename="{}.ic".format(gwtname) + ) + + # Instantiating MODFLOW 6 transport advection package + flopy.mf6.ModflowGwtadv( + gwt, scheme=scheme, pname="ADV-2", filename="{}.adv".format(gwtname) + ) + + # Instantiating MODFLOW 6 transport dispersion package + flopy.mf6.ModflowGwtdsp( + gwt, + xt3d_off=True, + alh=dispersivity, + alv=dispersivity, + ath1=dispersivity, + atv=dispersivity, + pname="CND-1", + filename="{}.cnd".format(gwtname), + ) + + # Instantiating MODFLOW 6 transport mass storage package (formerly "reaction" package in MT3DMS) + flopy.mf6.ModflowGwtmst( + gwt, + save_flows=True, + porosity=prsity, + sorption="linear", + bulk_density=rhob, + distcoef=Kd, + pname="MST-1", + filename="{}.mst".format(gwtname), + ) + + flopy.mf6.modflow.ModflowGwtuzt( + gwt, + boundnames=False, + save_flows=True, + print_input=True, + print_flows=True, + print_concentration=True, + concentration_filerecord=gwtname + ".uzt.bin", + budget_filerecord=gwtname + ".uzt.bud", + packagedata=uzt_pkdat, + uztperioddata=uzt_perdat, + flow_package_name="UZF-1", + pname="UZT-1", + filename="{}.uzt".format(gwtname), + ) + + flopy.mf6.modflow.ModflowGwtsft( + gwt, + boundnames=False, + flow_package_name="SFR-1", + print_concentration=True, + save_flows=True, + print_flows=True, + concentration_filerecord=gwtname + ".sft.bin", + budget_filerecord=gwtname + ".sft.bud", + packagedata=sft_pkdat, + reachperioddata=sft_perdat, + pname="SFT-1", + filename="{}.sft".format(gwtname), + ) + + flopy.mf6.modflow.ModflowGwtmvt( + gwt, + save_flows=True, + budget_filerecord=gwtname + ".mvt.bud", + filename="{}.mvt".format(gwtname), + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGwtoc( + gwt, + pname="OC-2", + budget_filerecord="{}.cbc".format(gwtname), + concentration_filerecord="{}.ucn".format(gwtname), + concentrationprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + ) + + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm( + gwt, sources=sourcerecarray, filename="{}.ssm".format(gwtname) + ) + + # Instantiating MODFLOW 6 flow-transport exchange mechanism + flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwtname, + pname="GWFGWT", + filename="{}.gwfgwt".format(gwtname), + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + + # read transport results from GWE model + name = cases[idx] + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + # overland flow + fpth = os.path.join(test.workspace, f"{gwfname}.mvr.bud") + + try: + # load mover flows + fobj = flopy.utils.CellBudgetFile(fpth, precision="double") + txtnames = fobj.get_unique_record_names() + times = fobj.get_times() + mvrdat = {} + for tm in times: + for txtname in txtnames: + # "MOVER-FLOW" should be the only text name + flowX = fobj.get_data(totim=tm, text=txtname) + mvrdat.update({tm: flowX}) + except: + assert False, f'could not load flow data from "{fpth}"' + + fpth = os.path.join(test.workspace, f"{gwtname}.mvt.bud") + try: + # load transport mover results + mobj = flopy.utils.CellBudgetFile(fpth, precision="double") + txtnames = mobj.get_unique_record_names() + times = fobj.get_times() + mvtdat = {} + for tm in times: + for txtname in txtnames: + massX = mobj.get_data(totim=tm, text=txtname) + mvtdat.update({tm: massX}) + except: + assert False, f'could not load transport data from "{fpth}"' + + # Mover amounts are representative of rejected infiltration only + # (there are no other MVR connections established). A known amount of + # water will be transferred based on the infiltration and VKS: 1 unit + # of water in stress period 1 and 2 units in the 2nd stress period + msg0 = "Rejected infiltration transfer amount not as expected" + for x in np.arange(len(mvrdat)): + for y in np.arange(len(mvrdat[x + 1])): + if len(mvrdat[x + 1][y]) == 0: + continue + else: + for z in np.arange(len(mvrdat[x + 1][y])): + assert np.isclose( + abs(mvrdat[x + 1][y][z][-1]), x + 1.0 + ), msg0 + + # Transport mover (MVT) amounts are known quantities + msg1 = "Rejected infiltration transfer mass amount not as expected" + for x in np.arange(len(mvtdat)): + for y in np.arange(len(mvtdat[x + 1])): + if len(mvtdat[x + 1][y]) == 0: + continue + else: + for z in np.arange(len(mvtdat[x + 1][y])): + assert np.isclose( + mvtdat[x + 1][y][z][-1], (x + 1.0) * concCell[z] + ), msg1 + + # Stream mass flows should be known + fpth = os.path.join(test.workspace, f"{gwtname}.sft.bud") + try: + # load mover flows + sfrobj = flopy.utils.CellBudgetFile(fpth, precision="double") + txtnames = sfrobj.get_unique_record_names() + txtnames = [nm.decode("utf-8").strip() for nm in txtnames] + times = sfrobj.get_times() + JAdat = {} + fromMvrDat = {} + for tm in times: + for txtname in txtnames: + if txtname == "FLOW-JA-FACE": + dat = sfrobj.get_data(totim=tm, text=txtname) + JAdat.update({tm: dat}) + if txtname == "FROM-MVR": + dat = sfrobj.get_data(totim=tm, text=txtname) + fromMvrDat.update({tm: dat}) + + except: + assert False, f'could not load flow data from "{fpth}"' + + msg2 = "Mass received by SFR ('FROM-MVR') not as expected" + for x in np.arange(len(fromMvrDat)): + for y in np.arange(len(fromMvrDat[x + 1][0])): + if fromMvrDat[x + 1][0][y][-1] == concCell[y]: + continue + else: + for z in np.arange(len(mvtdat[x + 1][y])): + assert np.isclose( + mvtdat[x + 1][y][z][-1], (x + 1.0) * concCell[z] + ), msg2 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(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/src/Utilities/BudgetTerm.f90 b/src/Utilities/BudgetTerm.f90 index 24de025e39d..4dd332dc89f 100644 --- a/src/Utilities/BudgetTerm.f90 +++ b/src/Utilities/BudgetTerm.f90 @@ -110,11 +110,13 @@ subroutine deallocate_arrays(this) ! -- dummy class(BudgetTermType) :: this ! - deallocate (this%id1) - deallocate (this%id2) - deallocate (this%flow) - deallocate (this%auxvar) - deallocate (this%auxtxt) + if (size(this%id1) > 0) then + deallocate (this%id1) + deallocate (this%id2) + deallocate (this%flow) + deallocate (this%auxvar) + deallocate (this%auxtxt) + end if end subroutine deallocate_arrays !> @brief reset the budget term and counter so terms can be updated From 402b64831e319b1584d8d2747f8a19a22a58181e Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Wed, 22 May 2024 13:37:02 -0700 Subject: [PATCH 2/8] trivial --- src/Utilities/BudgetObject.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Utilities/BudgetObject.f90 b/src/Utilities/BudgetObject.f90 index f3b03a3cbb0..0aba86b73de 100644 --- a/src/Utilities/BudgetObject.f90 +++ b/src/Utilities/BudgetObject.f90 @@ -482,7 +482,7 @@ subroutine save_flows(this, dis, ibinun, kstp, kper, delt, & real(DP), intent(in) :: pertim real(DP), intent(in) :: totim integer(I4B), intent(in) :: iout - ! -- dummy + ! -- local integer(I4B) :: i ! ! -- Save flows for each budget term From b493c9b43539f1615f570b40f1f72e6dd37311ac Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Wed, 22 May 2024 13:37:56 -0700 Subject: [PATCH 3/8] skip writing of flow-ja-face when only one reach --- src/Utilities/BudgetTerm.f90 | 38 +++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/Utilities/BudgetTerm.f90 b/src/Utilities/BudgetTerm.f90 index 4dd332dc89f..2722d49abe9 100644 --- a/src/Utilities/BudgetTerm.f90 +++ b/src/Utilities/BudgetTerm.f90 @@ -205,24 +205,26 @@ subroutine save_flows(this, dis, ibinun, kstp, kper, delt, pertim, totim, & end do ! ! -- Write the header - call ubdsv06(kstp, kper, this%flowtype, & - this%text1id1, this%text2id1, & - this%text1id2, this%text2id2, & - ibinun, this%naux, this%auxtxt, & - nlist, 1, 1, nlist, & - iout, delt, pertim, totim) - ! - ! -- Write each entry - do i = 1, this%nlist - q = this%flow(i) - n1 = this%id1(i) - n2 = this%id2(i) - if (n1 <= 0 .or. n2 <= 0) cycle - call dis%record_mf6_list_entry(ibinun, n1, n2, q, & - this%naux, this%auxvar(:, i), & - olconv=this%olconv1, & - olconv2=this%olconv2) - end do + if (nlist > 0) then + call ubdsv06(kstp, kper, this%flowtype, & + this%text1id1, this%text2id1, & + this%text1id2, this%text2id2, & + ibinun, this%naux, this%auxtxt, & + nlist, 1, 1, nlist, & + iout, delt, pertim, totim) + ! + ! -- Write each entry + do i = 1, this%nlist + q = this%flow(i) + n1 = this%id1(i) + n2 = this%id2(i) + if (n1 <= 0 .or. n2 <= 0) cycle + call dis%record_mf6_list_entry(ibinun, n1, n2, q, & + this%naux, this%auxvar(:, i), & + olconv=this%olconv1, & + olconv2=this%olconv2) + end do + end if ! end subroutine save_flows From 120cd67d5923d2a7217fe79cce92097c3aa7185f Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Wed, 22 May 2024 16:16:09 -0700 Subject: [PATCH 4/8] stop writing blank entries to output, patch up some autotests as a result --- autotest/test_gwf_mvr01.py | 55 ++++-------- autotest/test_gwt_uztmvt2x2.py | 13 ++- src/Model/GroundWaterFlow/gwf-sfr.f90 | 120 ++++++++++++++------------ 3 files changed, 87 insertions(+), 101 deletions(-) diff --git a/autotest/test_gwf_mvr01.py b/autotest/test_gwf_mvr01.py index 08368bf553f..a6ab48a20a8 100644 --- a/autotest/test_gwf_mvr01.py +++ b/autotest/test_gwf_mvr01.py @@ -374,48 +374,29 @@ def check_output(idx, test): times = bobj.get_times() records = bobj.get_data(totim=times[-1]) adt = [("node", " @ brief Output package flow terms. !! - !! Output SFR package flow terms. - !! + !! Write flows to binary file and/or print flows to budget !< subroutine sfr_ot_package_flows(this, icbcfl, ibudfl) ! -- modules @@ -2623,7 +2622,7 @@ subroutine sfr_ot_package_flows(this, icbcfl, ibudfl) ! -- dummy variables class(SfrType) :: this !< SfrType object integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output - integer(I4B), intent(in) :: ibudfl !< flag indication if cell-by-cell data should be saved + integer(I4B), intent(in) :: ibudfl !< flag indicating if cell-by-cell data should be saved ! -- local variables integer(I4B) :: ibinun character(len=20), dimension(:), allocatable :: cellidstr @@ -5097,6 +5096,11 @@ subroutine sfr_setup_budobj(this) ! so they can be written to the binary budget files, but these internal ! flows are not included as part of the budget table. nbudterm = 8 + ! + ! -- GWF models with a single reach should not have an entry + if (this%nconn == 0) nbudterm = nbudterm - 1 + ! + ! -- Account for mover connection and aux vars in budget object if (this%imover == 1) nbudterm = nbudterm + 2 if (this%naux > 0) nbudterm = nbudterm + 1 ! @@ -5107,29 +5111,33 @@ subroutine sfr_setup_budobj(this) idx = 0 ! ! -- Go through and set up each budget term - text = ' FLOW-JA-FACE' - idx = idx + 1 - maxlist = this%nconn - naux = 1 - auxtxt(1) = ' FLOW-AREA' - call this%budobj%budterm(idx)%initialize(text, & - this%name_model, & - this%packName, & - this%name_model, & - this%packName, & - maxlist, .false., .false., & - naux, auxtxt) ! - ! -- store connectivity - call this%budobj%budterm(idx)%reset(this%nconn) - q = DZERO - do n = 1, this%maxbound - n1 = n - do i = this%ia(n) + 1, this%ia(n + 1) - 1 - n2 = this%ja(i) - call this%budobj%budterm(idx)%update_term(n1, n2, q) + ! -- GWF models with a single reach do not have flow-ja-face connections + if (this%nconn /= 0) then + text = ' FLOW-JA-FACE' + idx = idx + 1 + maxlist = this%nconn + naux = 1 + auxtxt(1) = ' FLOW-AREA' + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux, auxtxt) + ! + ! -- store connectivity + call this%budobj%budterm(idx)%reset(this%nconn) + q = DZERO + do n = 1, this%maxbound + n1 = n + do i = this%ia(n) + 1, this%ia(n + 1) - 1 + n2 = this%ja(i) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + end do end do - end do + end if ! ! -- text = ' GWF' @@ -5317,40 +5325,42 @@ subroutine sfr_fill_budobj(this) idx = 0 ! ! -- FLOW JA FACE - idx = idx + 1 - call this%budobj%budterm(idx)%reset(this%nconn) - do n = 1, this%maxbound - n1 = n - q = DZERO - ca = DZERO - do i = this%ia(n) + 1, this%ia(n + 1) - 1 - n2 = this%ja(i) - if (this%iboundpak(n) /= 0) then - ! flow to downstream reaches - if (this%idir(i) < 0) then - qt = this%dsflow(n) - q = -this%qconn(i) - ! flow from upstream reaches + if (this%nconn /= 0) then + idx = idx + 1 + call this%budobj%budterm(idx)%reset(this%nconn) + do n = 1, this%maxbound + n1 = n + q = DZERO + ca = DZERO + do i = this%ia(n) + 1, this%ia(n + 1) - 1 + n2 = this%ja(i) + if (this%iboundpak(n) /= 0) then + ! flow to downstream reaches + if (this%idir(i) < 0) then + qt = this%dsflow(n) + q = -this%qconn(i) + ! flow from upstream reaches + else + qt = this%usflow(n) + do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1 + if (this%idir(ii) > 0) cycle + if (this%ja(ii) /= n) cycle + q = this%qconn(ii) + exit + end do + end if + ! calculate flow area + call this%sfr_calc_reach_depth(n, qt, d) + ca = this%calc_area_wet(n, d) else - qt = this%usflow(n) - do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1 - if (this%idir(ii) > 0) cycle - if (this%ja(ii) /= n) cycle - q = this%qconn(ii) - exit - end do + q = DZERO + ca = DZERO end if - ! calculate flow area - call this%sfr_calc_reach_depth(n, qt, d) - ca = this%calc_area_wet(n, d) - else - q = DZERO - ca = DZERO - end if - this%qauxcbc(1) = ca - call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc) + this%qauxcbc(1) = ca + call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc) + end do end do - end do + end if ! ! -- GWF (LEAKAGE) idx = idx + 1 From 73da87c9e93e18aff192aea10e2798addbd67fe8 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 8 Jul 2024 07:43:59 -0700 Subject: [PATCH 5/8] update a new autotest to keep inline with refactor of gweest input --- autotest/test_gwe_drycell_cnd3.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/autotest/test_gwe_drycell_cnd3.py b/autotest/test_gwe_drycell_cnd3.py index 0a72cedb324..67d0b922145 100644 --- a/autotest/test_gwe_drycell_cnd3.py +++ b/autotest/test_gwe_drycell_cnd3.py @@ -341,9 +341,11 @@ def add_gwe_model(sim, gwename): gwe, save_flows=True, porosity=prsity, + heat_capacity_water=cpw, + density_water=rhow, + latent_heat_vaporization=lhv, cps=cps, rhos=rhos, - packagedata=[cpw, rhow, lhv], pname="EST", filename="{}.est".format(gwename), ) From 113f17448a89103a3e68cf72b1a8066a16fe79f4 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 8 Jul 2024 09:56:10 -0700 Subject: [PATCH 6/8] some regression comparison tests not expected to pass with this PR --- autotest/test_testmodels_mf6.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/autotest/test_testmodels_mf6.py b/autotest/test_testmodels_mf6.py index 52a923d2f90..bb585ce60ae 100644 --- a/autotest/test_testmodels_mf6.py +++ b/autotest/test_testmodels_mf6.py @@ -26,6 +26,11 @@ "test014_NWTP3Low_dev", ] +excluded_comparisons = [ + "test051_uzfp3_lakmvr_v2_dev", + "test051_uzfp3_wellakmvr_v2", +] + @pytest.mark.repo @pytest.mark.regression @@ -39,13 +44,18 @@ def test_model( ): model_path = test_model_mf6.parent model_name = model_path.name - excluded = model_name in excluded_models + exclude_model = model_name in excluded_models + exclude_comparison = model_name in excluded_comparisons compare = ( - get_mf6_comparison(model_path) if original_regression else "mf6_regression" + get_mf6_comparison(model_path) + if original_regression + else None + if exclude_comparison + else "mf6_regression" ) dev_only = "dev" in model_name and "not developmode" in markers - if excluded or dev_only: - reason = "excluded" if excluded else "developmode only" + if exclude_model or dev_only: + reason = "excluded" if exclude_model else "developmode only" pytest.skip(f"Skipping: {model_name} ({reason})") # setup test workspace and framework From b9b8a9760e9ed5175be730aaf17a08f72110095b Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 30 Sep 2024 12:17:15 -0700 Subject: [PATCH 7/8] get up to date with ruff --- autotest/test_gwe_drycell_cnd3.py | 9 +++------ autotest/test_gwf_sfrsft_gwdischrg.py | 8 ++++---- autotest/test_gwt_uztmvt2x1.py | 7 +++---- 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/autotest/test_gwe_drycell_cnd3.py b/autotest/test_gwe_drycell_cnd3.py index 67d0b922145..36f88996ba1 100644 --- a/autotest/test_gwe_drycell_cnd3.py +++ b/autotest/test_gwe_drycell_cnd3.py @@ -25,10 +25,10 @@ # Imports import os + +import flopy import numpy as np import pytest -import flopy - from framework import TestFramework @@ -145,7 +145,6 @@ def isMonotonic(A): def add_gwf_model(sim, gwfname, newton=False): - # Instantiating MODFLOW 6 groundwater flow model if newton: gwf = flopy.mf6.ModflowGwf( @@ -275,7 +274,6 @@ def add_gwf_model(sim, gwfname, newton=False): def add_gwe_model(sim, gwename): - gwe = flopy.mf6.ModflowGwe( sim, modelname=gwename, model_nam_file="{}.nam".format(gwename) ) @@ -402,12 +400,11 @@ def add_gwe_model(sim, gwename): def build_models(idx, test): - # Base MF6 GWF model type ws = test.workspace name = cases[idx] - print("Building MF6 model...()".format(name)) + print("Building MF6 model...{}".format(name)) # generate names for each model gwfname1 = "gwf-" + name + "nwt1" diff --git a/autotest/test_gwf_sfrsft_gwdischrg.py b/autotest/test_gwf_sfrsft_gwdischrg.py index d905697b292..11988a3d2e8 100644 --- a/autotest/test_gwf_sfrsft_gwdischrg.py +++ b/autotest/test_gwf_sfrsft_gwdischrg.py @@ -1,13 +1,11 @@ # Test groundwater discharge to a stream and then go on to test # that a transport model with a single reach works. -import math import pathlib as pl import flopy import numpy as np import pytest - from framework import TestFramework cases = ["sfr-gwfout", "sfr-gwf-trnsprt"] @@ -224,7 +222,8 @@ def check_output(idx, test): try: # load simulated concentration in SFT cobj = flopy.utils.HeadFile( - sft_obs_fl, text="CONCENTRATION" # precision="double" + sft_obs_fl, + text="CONCENTRATION", # precision="double" ) sim_conc_sft = cobj.get_alldata() except: @@ -236,7 +235,8 @@ def check_output(idx, test): try: # load simulated concentration of groundwater cobj = flopy.utils.HeadFile( - gwt_sim_conc, text="CONCENTRATION" # precision="double" + gwt_sim_conc, + text="CONCENTRATION", # precision="double" ) conc_gw = cobj.get_alldata() except: diff --git a/autotest/test_gwt_uztmvt2x1.py b/autotest/test_gwt_uztmvt2x1.py index 0b78660ebfa..c0b468344bc 100644 --- a/autotest/test_gwt_uztmvt2x1.py +++ b/autotest/test_gwt_uztmvt2x1.py @@ -33,10 +33,10 @@ # Imports import os + +import flopy import numpy as np import pytest -import flopy - from framework import TestFramework # Base simulation and model name and workspace @@ -177,12 +177,11 @@ def build_models(idx, test): - # Base MF6 GWF model type ws = test.workspace name = cases[idx] - print("Building MF6 model...()".format(name)) + print("Building MF6 model...{}".format(name)) # generate names for each model gwfname = "gwf-" + name From b28c6319cda07ba8742b890747ecc37fd0f1984f Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Tue, 5 Nov 2024 13:24:14 -0800 Subject: [PATCH 8/8] update argument usage --- autotest/test_gwe_drycell_cnd3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autotest/test_gwe_drycell_cnd3.py b/autotest/test_gwe_drycell_cnd3.py index 36f88996ba1..dc1515c2cb6 100644 --- a/autotest/test_gwe_drycell_cnd3.py +++ b/autotest/test_gwe_drycell_cnd3.py @@ -342,8 +342,8 @@ def add_gwe_model(sim, gwename): heat_capacity_water=cpw, density_water=rhow, latent_heat_vaporization=lhv, - cps=cps, - rhos=rhos, + heat_capacity_solid=cps, + density_solid=rhos, pname="EST", filename="{}.est".format(gwename), )