diff --git a/autotest/test_netcdf_gwe_cnd.py b/autotest/test_netcdf_gwe_cnd.py index 46589e2ca04..02a542f31af 100644 --- a/autotest/test_netcdf_gwe_cnd.py +++ b/autotest/test_netcdf_gwe_cnd.py @@ -1,14 +1,7 @@ """ -Test problem for GWE - -One-Dimensional Transport in a Uniform Flow Field. -The purpose of this script is to test the new heat transport model developed -for MODFLOW 6. To that end, this problem uses the setup of the first MT3DMS -test problem but adapts it for heat. MODFLOW 6 is setup using the new GWE -model with input parameters entered in their native units. - -It may be possible to find a 1D heat transport analytical solution in the -future. +NetCDF export test version of test_gwe_cnd. This test compares +the temperature and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. """ # Imports @@ -17,20 +10,6 @@ import numpy as np import pytest -try: - import xarray as xa - import xugrid as xu -except ImportError: - pytest.skip("xuarray and xugrid not found", allow_module_level=True) - -try: - import pymake -except: - msg = "Error. Pymake package is not available.\n" - msg += "Try installing using the following command:\n" - msg += " pip install https://github.com/modflowpy/pymake/zipball/master" - raise Exception(msg) - try: import flopy except: @@ -39,332 +18,34 @@ msg += " pip install flopy" raise Exception(msg) +try: + import xarray as xa + import xugrid as xu +except ImportError: + pytest.skip("xuarray and xugrid not found", allow_module_level=True) from framework import TestFramework -# Base simulation and model name and workspace - -viscosity_on = [False] cases = ["cnd01"] -# Model units - -length_units = "meters" -time_units = "days" - -# Table MODFLOW 6 GWE comparison to MT3DMS - -nper = 1 # Number of periods -nlay = 1 # Number of layers -ncol = 101 # Number of columns -nrow = 1 # Number of rows -delr = 10.0 # Column width ($m$) -delc = 1.0 # Row width ($m$) -top = 0.0 # Top of the model ($m$) -botm = -1.0 # Layer bottom elevations ($m$) -prsity = 0.25 # Porosity -perlen = 2000 # Simulation time ($days$) -k11 = 1.0 # Horizontal hydraulic conductivity ($m/d$) - -# Set some static model parameter values - -k33 = k11 # Vertical hydraulic conductivity ($m/d$) -laytyp = 1 -nstp = 100.0 -dt0 = perlen / nstp -Lx = (ncol - 1) * delr -v = 0.24 -q = v * prsity -h1 = q * Lx -strt = np.zeros((nlay, nrow, ncol), dtype=float) -strt[0, 0, 0] = h1 # Starting head ($m$) -l = 1000.0 # Needed for plots -icelltype = 1 # Cell conversion type -ibound = np.ones((nlay, nrow, ncol), dtype=int) -ibound[0, 0, 0] = -1 -ibound[0, 0, -1] = -1 - -# Set some static transport related model parameter values - -mixelm = 0 # FD -rhob = 1110.0 -sp2 = 0.0 # read, but not used in this problem -kd = 1.8168e-4 -strt_temp = np.zeros((nlay, nrow, ncol), dtype=float) -dispersivity = 1.0 -dmcoef = 3.2519e-7 # Molecular diffusion coefficient - -# Set some static heat transport related model parameter values -cpw = 4183.0 -rhow = 1000.0 -lhv = 2454.0 - -# Set solver parameter values (and related) -nouter, ninner = 100, 300 -hclose, rclose, relax = 1e-6, 1e-6, 1.0 -ttsmult = 1.0 -dceps = 1.0e-5 # HMOC parameters in case they are invoked -nplane = 1 # HMOC -npl = 0 # HMOC -nph = 4 # HMOC -npmin = 0 # HMOC -npmax = 8 # HMOC -nlsink = nplane # HMOC -npsink = nph # HMOC - -# Static temporal data used by TDIS file - -tdis_rc = [] -tdis_rc.append((perlen, nstp, 1.0)) - -# ### Create MODFLOW 6 GWE MT3DMS Example 1 Boundary Conditions -# -# Constant head cells are specified on both ends of the model - -chdspd = [[(0, 0, 0), h1], [(0, 0, ncol - 1), 0.0]] -c0 = 40.0 -ctpspd = [[(0, 0, 0), c0]] - def build_models(idx, test): - # Base MF6 GWE model type - ws = test.workspace - name = cases[idx] - - print("Building MF6 model...()".format(name)) - - # generate names for each model - gwfname = "gwf-" + name - gwename = "gwe-" + name - - sim_ws = os.path.join(ws, 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, start_date_time="2041-01-01T00:00:00-05:00", - ) + from test_gwe_cnd import build_models as build - # 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=botm, - idomain=np.ones((nlay, nrow, ncol), dtype=int), - filename="{}.dis".format(gwfname), - ) - - # Instantiating MODFLOW 6 node-property flow package - flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=False, - icelltype=icelltype, - k=k11, - k33=k33, - save_specific_discharge=True, - filename="{}.npf".format(gwfname), - ) - - # Instantiating MODFLOW 6 initial conditions package for flow model - flopy.mf6.ModflowGwfic(gwf, strt=strt, filename="{}.ic".format(gwfname)) - - # Instantiating VSC - if viscosity_on[idx]: - # Instantiate viscosity (VSC) package - vsc_filerecord = "{}.vsc.bin".format(gwfname) - vsc_pd = [(0, 0.0, 20.0, gwename, "temperature")] - flopy.mf6.ModflowGwfvsc( - gwf, - viscref=8.904e-4, - viscosity_filerecord=vsc_filerecord, - thermal_formulation="nonlinear", - thermal_a2=10.0, - thermal_a3=248.37, - thermal_a4=133.16, - nviscspecies=len(vsc_pd), - packagedata=vsc_pd, - pname="vsc", - filename="{}.vsc".format(gwfname), - ) - - # Instantiating MODFLOW 6 constant head package - flopy.mf6.ModflowGwfchd( - gwf, - maxbound=len(chdspd), - stress_period_data=chdspd, - save_flows=False, - pname="CHD-1", - filename="{}.chd".format(gwfname), - ) - - # Instantiating MODFLOW 6 output control package for flow model - flopy.mf6.ModflowGwfoc( - gwf, - 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 groundwater transport package - gwe = flopy.mf6.MFModel( - sim, - model_type="gwe6", - modelname=gwename, - model_nam_file="{}.nam".format(gwename), - ) - gwe.name_file.save_flows = True + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwe = sim.gwe[0] gwe.name_file.export_netcdf = "ugrid" - 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, - export_array_netcdf=True, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - idomain=1, - filename="{}.dis".format(gwename), - ) - - # Instantiating MODFLOW 6 transport initial concentrations - flopy.mf6.ModflowGweic( - gwe, export_array_netcdf=True, strt=strt_temp, filename="{}.ic".format(gwename) - ) - - # Instantiating MODFLOW 6 transport advection package - if mixelm == 0: - scheme = "UPSTREAM" - elif mixelm == -1: - scheme = "TVD" - else: - raise Exception() - flopy.mf6.ModflowGweadv( - gwe, scheme=scheme, filename="{}.adv".format(gwename) - ) - - # Instantiating MODFLOW 6 transport dispersion package - if dispersivity != 0: - flopy.mf6.ModflowGwecnd( - gwe, - xt3d_off=True, - export_array_netcdf=True, - alh=dispersivity, - ath1=dispersivity, - ktw=0.5918, - kts=0.2700, - filename="{}.cnd".format(gwename), - ) - - # Instantiating MODFLOW 6 transport mass storage package (formerly "reaction" package in MT3DMS) - flopy.mf6.ModflowGweest( - gwe, - save_flows=True, - heat_capacity_water=cpw, - density_water=rhow, - porosity=prsity, - cps=760.0, - rhos=1500.0, - filename="{}.est".format(gwename), - ) - - # Instantiating MODFLOW 6 transport constant concentration package - flopy.mf6.ModflowGwectp( - gwe, - maxbound=len(ctpspd), - stress_period_data=ctpspd, - save_flows=False, - pname="CTP-1", - filename="{}.ctp".format(gwename), - ) - - # Instantiating MODFLOW 6 transport source-sink mixing package - flopy.mf6.ModflowGwessm( - gwe, sources=[[]], filename="{}.ssm".format(gwename) - ) - - # Instantiate MODFLOW 6 heat transport output control package - flopy.mf6.ModflowGweoc( - gwe, - 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")], - ) - - # Instantiating MODFLOW 6 flow-transport exchange mechanism - flopy.mf6.ModflowGwfgwe( - sim, - exgtype="GWF6-GWE6", - exgmnamea=gwfname, - exgmnameb=gwename, - filename="{}.gwfgwe".format(name), - ) - - return sim, None + gwe.dis.export_array_netcdf = True + gwe.ic.export_array_netcdf = True + gwe.cnd.export_array_netcdf = True + return sim, dummy def check_output(idx, test): - print("evaluating results...") + from test_gwe_cnd import check_output as check + + check(idx, test) # read transport results from GWE model name = cases[idx] @@ -380,114 +61,7 @@ def check_output(idx, test): except: assert False, f'could not load concentration data from "{fpth}"' - # This is the answer - c_ans = [ - 4.00000000e01, - 3.99999983e01, - 3.99999898e01, - 3.99999566e01, - 3.99998462e01, - 3.99995197e01, - 3.99986427e01, - 3.99964775e01, - 3.99915230e01, - 3.99809477e01, - 3.99597839e01, - 3.99198995e01, - 3.98488519e01, - 3.97288247e01, - 3.95359427e01, - 3.92403042e01, - 3.88070317e01, - 3.81985089e01, - 3.73777505e01, - 3.63125911e01, - 3.49801399e01, - 3.33708033e01, - 3.14911723e01, - 2.93652158e01, - 2.70334931e01, - 2.45504338e01, - 2.19800532e01, - 1.93907148e01, - 1.68496655e01, - 1.44180473e01, - 1.21469471e01, - 1.00748333e01, - 8.22648357e00, - 6.61329449e00, - 5.23470060e00, - 4.08034410e00, - 3.13261741e00, - 2.36924164e00, - 1.76562010e00, - 1.29679741e00, - 9.38944408e-01, - 6.70362685e-01, - 4.72056032e-01, - 3.27947150e-01, - 2.24829892e-01, - 1.52144844e-01, - 1.01654320e-01, - 6.70766201e-02, - 4.37223104e-02, - 2.81598160e-02, - 1.79249349e-02, - 1.12795213e-02, - 7.01828727e-03, - 4.31895689e-03, - 2.62924728e-03, - 1.58374083e-03, - 9.44125798e-04, - 5.57133590e-04, - 3.25507431e-04, - 1.88330495e-04, - 1.07925092e-04, - 6.12700035e-05, - 3.44648666e-05, - 1.92125906e-05, - 1.06157638e-05, - 5.81494908e-06, - 3.15821246e-06, - 1.70101068e-06, - 9.08679391e-07, - 4.81524218e-07, - 2.53159103e-07, - 1.32068539e-07, - 6.83748562e-08, - 3.51353218e-08, - 1.79225415e-08, - 9.07652498e-09, - 4.56413759e-09, - 2.27913640e-09, - 1.13033292e-09, - 5.56823550e-10, - 2.72491770e-10, - 1.32483548e-10, - 6.40015158e-11, - 3.07244529e-11, - 1.46584603e-11, - 6.95098705e-12, - 3.27643160e-12, - 1.53530190e-12, - 7.15261898e-13, - 3.31325318e-13, - 1.52616350e-13, - 6.99104644e-14, - 3.18504005e-14, - 1.44329547e-14, - 6.50576657e-15, - 2.91728603e-15, - 1.30145909e-15, - 5.77678170e-16, - 2.55141072e-16, - 1.12178999e-16, - 5.01900830e-17, - ] - - msg = f"gwe temperatures do not match stored concentrations" - assert np.allclose(conc1[-1, 0, 0, :], c_ans, atol=1e-5), msg - + # Check NetCDF output nc_fpth = os.path.join(test.workspace, f"{gwename}.nc") ds = xu.open_dataset(nc_fpth) xds = ds.ugrid.to_dataset() @@ -552,7 +126,7 @@ def check_output(idx, test): ), f"NetCDF input array comparison failure, variable={var}" -# - No need to change any code below +@pytest.mark.netcdf @pytest.mark.parametrize( "idx, name", list(enumerate(cases)), diff --git a/autotest/test_netcdf_gwf_disv.py b/autotest/test_netcdf_gwf_disv.py index f44cc546644..a65634ca8fe 100644 --- a/autotest/test_netcdf_gwf_disv.py +++ b/autotest/test_netcdf_gwf_disv.py @@ -1,10 +1,7 @@ """ -Test of GWF DISV Package. Use the flopy disv tool to create -a simple regular grid example, but using DISV instead of DIS. -Use a large offset for x and y vertices to ensure that the area -calculation in MODFLOW 6 is correct. For the second case, set -one of the cells inactive and test to make sure connectivity -in binary grid file is correct. +NetCDF export test version of test_gwf_disv. This test compares +the heads and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. """ import os @@ -22,14 +19,14 @@ except ImportError: pytest.skip("xuarray and xugrid not found", allow_module_level=True) -cases = ["disv01a_ncf", "disv01b_ncf"] +cases = ["disv01a", "disv01b"] wkt = ( 'PROJCS["NAD83 / UTM zone 18N", ' 'GEOGCS["NAD83", ' 'DATUM["North_American_Datum_1983", ' 'SPHEROID["GRS 1980",6378137,298.257222101], ' - 'TOWGS84[0,0,0,0,0,0,0]], ' + "TOWGS84[0,0,0,0,0,0,0]], " 'PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], ' 'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], ' 'AUTHORITY["EPSG","4269"]], ' @@ -47,86 +44,40 @@ def build_models(idx, test): + from test_gwf_disv import build_models as build + + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = "ugrid" + gwf.disv.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True + name = cases[idx] - ws = test.workspace - nlay = 3 - nrow = 3 - ncol = 3 - delr = 10.0 - delc = 10.0 - top = 0 - botm = [-10, -20, -30] - xoff = 100000000.0 - yoff = 100000000.0 - disvkwargs = get_disv_kwargs( - nlay, - nrow, - ncol, - delr, - delc, - top, - botm, - xoff, - yoff, - ) - disvkwargs["export_array_netcdf"] = True - if idx == 1: - # for the second test, set one cell to idomain = 0 - idomain = np.ones((nlay, nrow * ncol), dtype=int) - idomain[0, 1] = 0 - disvkwargs["idomain"] = idomain - - sim = flopy.mf6.MFSimulation( - sim_name=name, - version="mf6", - exe_name="mf6", - sim_ws=ws, - ) - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", start_date_time="2041-01-01T00:00:00-05:00" + + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.disv, ogc_wkt=wkt, filename=f"{name}.disv.ncf" ) - gwf = flopy.mf6.ModflowGwf(sim, modelname=name, export_netcdf="ugrid") - ims = flopy.mf6.ModflowIms(sim, print_option="SUMMARY") - disv = flopy.mf6.ModflowGwfdisv(gwf, **disvkwargs) - # netcdf configuration - ncf = flopy.mf6.ModflowUtlncf(disv, ogc_wkt=wkt, filename=f"{name}.disv.ncf") - ic = flopy.mf6.ModflowGwfic(gwf, export_array_netcdf=True, strt=0.0) - npf = flopy.mf6.ModflowGwfnpf(gwf, export_array_netcdf=True) - spd = {0: [[(0, 0), 1.0], [(0, nrow * ncol - 1), 0.0]]} - chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd) # output control oc = flopy.mf6.ModflowGwfoc( gwf, - # budget_filerecord=f"{name}.cbc", head_filerecord=f"{name}.hds", - # headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], - # saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], saverecord=[ ("HEAD", "ALL"), ], - # printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")], ) - return sim, None + return sim, dummy def check_output(idx, test): - name = test.name + from test_gwf_disv import check_output as check - fname = os.path.join(test.workspace, name + ".disv.grb") - grbobj = flopy.mf6.utils.MfGrdFile(fname) - ncpl = grbobj._datadict["NCPL"] - ia = grbobj._datadict["IA"] - ja = grbobj._datadict["JA"] - - if idx == 1: - # assert ncpl == disvkwargs["ncpl"] - assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7])) - assert np.array_equal(ja[:6], np.array([1, 4, 10, 3, 6, 12])) - assert ia[-1] == 127 - assert ia.shape[0] == 28, "ia should have size of 28" - assert ja.shape[0] == 126, "ja should have size of 126" + check(idx, test) # Check NetCDF output + name = test.name nc_fpth = os.path.join(test.workspace, name + ".nc") ds = xu.open_dataset(nc_fpth) xds = ds.ugrid.to_dataset() @@ -148,8 +99,9 @@ def check_output(idx, test): for l in range(nlay): assert np.allclose( np.array(rec[l]).flatten(), - #xds[f"head_l{l+1}"][timestep, :].data, - xds[f"head_l{l+1}"][timestep, :].fillna(1.00000000e+30).data, + xds[f"head_l{l+1}"][timestep, :] + .fillna(1.00000000e30) + .data, ), f"NetCDF-Headfile comparison failure in timestep {timestep+1}" timestep += 1 diff --git a/autotest/test_netcdf_gwf_sto01.py b/autotest/test_netcdf_gwf_sto01.py index fee90e12f4b..802bdff6c4a 100644 --- a/autotest/test_netcdf_gwf_sto01.py +++ b/autotest/test_netcdf_gwf_sto01.py @@ -1,3 +1,7 @@ +""" +NetCDF export test version of test_gwf_sto01. +""" + import os import flopy @@ -13,100 +17,15 @@ except ImportError: pytest.skip("xuarray and xugrid not found", allow_module_level=True) -cases = ["gwf_sto01_ncf"] -cmppth = "mfnwt" -tops = [0.0] +cases = ["gwf_sto01"] htol = [None for _ in range(len(cases))] -dtol = 1e-3 -budtol = 1e-2 -bud_lst = ["STO-SS_IN", "STO-SS_OUT", "STO-SY_IN", "STO-SY_OUT"] - -# static model data -# temporal discretization -nper = 31 -perlen = [1.0] + [365.2500000 for _ in range(nper - 1)] -nstp = [1] + [6 for _ in range(nper - 1)] -tsmult = [1.0] + [1.3 for _ in range(nper - 1)] -# tsmult = [1.0] + [1.0 for i in range(nper - 1)] -steady = [True] + [False for _ in range(nper - 1)] -tdis_rc = [] -for i in range(nper): - tdis_rc.append((perlen[i], nstp[i], tsmult[i])) - -# spatial discretization data -nlay, nrow, ncol = 3, 10, 10 -shape3d = (nlay, nrow, ncol) -size3d = nlay * nrow * ncol -delr, delc = 1000.0, 2000.0 -botm = [-100, -150.0, -350.0] -strt = 0.0 -hnoflo = 1e30 -hdry = -1e30 - -# calculate hk -hk1fact = 1.0 / 50.0 -hk1 = np.ones((nrow, ncol), dtype=float) * 0.5 * hk1fact -hk1[0, :] = 1000.0 * hk1fact -hk1[-1, :] = 1000.0 * hk1fact -hk1[:, 0] = 1000.0 * hk1fact -hk1[:, -1] = 1000.0 * hk1fact -hk = [20.0, hk1, 5.0] - -# calculate vka -vka = [1e6, 7.5e-5, 1e6] - -# all cells are active and layer 1 is convertible -ib = 1 -laytyp = [1, 0, 0] - -# solver options -nouter, ninner = 500, 300 -hclose, rclose, relax = 1e-9, 1e-6, 1.0 -newtonoptions = "NEWTON" -imsla = "BICGSTAB" - -# chd data -c = [] -c6 = [] -ccol = [3, 4, 5, 6] -for j in ccol: - c.append([0, nrow - 1, j, strt, strt]) - c6.append([(0, nrow - 1, j), strt]) -cd = {0: c} -cd6 = {0: c6} -maxchd = len(cd[0]) - -# pumping well data -wr = [0, 0, 0, 0, 1, 1, 2, 2, 3, 3] -wc = [0, 1, 8, 9, 0, 9, 0, 9, 0, 0] -wrp = [2, 2, 3, 3] -wcp = [5, 6, 5, 6] -wq = [-14000.0, -8000.0, -5000.0, -3000.0] -d = [] -d6 = [] -for r, c, q in zip(wrp, wcp, wq): - d.append([2, r, c, q]) - d6.append([(2, r, c), q]) -wd = {1: d} -wd6 = {1: d6} -maxwel = len(wd[1]) - -# recharge data -q = 3000.0 / (delr * delc) -v = np.zeros((nrow, ncol), dtype=float) -for r, c in zip(wr, wc): - v[r, c] = q -rech = {0: v} - -# storage and compaction data -ske = [6e-4, 3e-4, 6e-4] wkt = ( 'PROJCS["NAD83 / UTM zone 18N", ' 'GEOGCS["NAD83", ' 'DATUM["North_American_Datum_1983", ' 'SPHEROID["GRS 1980",6378137,298.257222101], ' - 'TOWGS84[0,0,0,0,0,0,0]], ' + "TOWGS84[0,0,0,0,0,0,0]], " 'PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], ' 'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], ' 'AUTHORITY["EPSG","4269"]], ' @@ -123,268 +42,30 @@ ) -# variant SUB package problem 3 def build_models(idx, test): - name = cases[idx] - - # build MODFLOW 6 files - ws = test.workspace - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws - ) - # create tdis package - tdis = flopy.mf6.ModflowTdis( - sim, - time_units="DAYS", - start_date_time="2041-01-01T00:00:00-05:00", - nper=nper, - perioddata=tdis_rc, - ) - - # create gwf model - top = tops[idx] - zthick = [top - botm[0], botm[0] - botm[1], botm[1] - botm[2]] - elevs = [top] + botm - - gwf = flopy.mf6.ModflowGwf( - sim, - modelname=name, - newtonoptions=newtonoptions, - save_flows=True, - export_netcdf="ugrid", - ) - - # create iterative model solution and register the gwf model with it - ims = 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=imsla, - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=relax, - ) - sim.register_ims_package(ims, [gwf.name]) - - dis = flopy.mf6.ModflowGwfdis( - gwf, - export_array_netcdf=True, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - filename=f"{name}.dis", - ) - - # netcdf configuration - flopy.mf6.ModflowUtlncf( - dis, - ogc_wkt=wkt, - filename=f"{name}.dis.ncf", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic( - gwf, export_array_netcdf=True, strt=strt, filename=f"{name}.ic" - ) - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=False, - export_array_netcdf=True, - # dev_modflowusg_upstream_weighted_saturation=True, - icelltype=laytyp, - k=hk, - k33=vka, - ) - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=False, - iconvert=laytyp, - ss=ske, - sy=0, - storagecoefficient=None, - steady_state={0: True}, - transient={1: True}, - ) + from test_gwf_sto01 import build_models as build - # recharge - rch = flopy.mf6.ModflowGwfrcha(gwf, readasarrays=True, recharge=rech) + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = "ugrid" + gwf.dis.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True - # wel file - wel = flopy.mf6.ModflowGwfwel( - gwf, - print_input=True, - print_flows=True, - maxbound=maxwel, - stress_period_data=wd6, - save_flows=False, - ) + name = cases[idx] - # chd files - chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd( - gwf, maxbound=maxchd, stress_period_data=cd6, save_flows=False + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.dis, ogc_wkt=wkt, filename=f"{name}.dis.ncf" ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.cbc", - head_filerecord=f"{name}.hds", - headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], - saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")], - ) - - # build MODFLOW-NWT files - cpth = cmppth - ws = os.path.join(test.workspace, cpth) - mc = flopy.modflow.Modflow( - name, - model_ws=ws, - version=cpth, - exe_name=try_get_target(test.targets, "mfnwt"), - ) - dis = flopy.modflow.ModflowDis( - mc, - nlay=nlay, - nrow=nrow, - ncol=ncol, - nper=nper, - perlen=perlen, - nstp=nstp, - tsmult=tsmult, - steady=steady, - delr=delr, - delc=delc, - top=top, - botm=botm, - ) - bas = flopy.modflow.ModflowBas( - mc, ibound=ib, strt=strt, hnoflo=hnoflo, stoper=0.01 - ) - upw = flopy.modflow.ModflowUpw( - mc, - laytyp=laytyp, - ipakcb=1001, - hk=hk, - vka=vka, - ss=ske, - sy=0.0, - hdry=hdry, - ) - chd = flopy.modflow.ModflowChd(mc, stress_period_data=cd) - rch = flopy.modflow.ModflowRch(mc, rech=rech) - wel = flopy.modflow.ModflowWel(mc, stress_period_data=wd) - oc = flopy.modflow.ModflowOc( - mc, - stress_period_data=None, - save_every=1, - save_types=["save head", "save budget"], - ) - fluxtol = (float(nlay * nrow * ncol) - 4.0) * rclose - nwt = flopy.modflow.ModflowNwt( - mc, - headtol=hclose, - fluxtol=fluxtol, - maxiterout=nouter, - linmeth=2, - maxitinner=ninner, - unitnumber=132, - options="SPECIFIED", - backflag=0, - idroptol=0, - ) - return sim, mc + return sim, dummy def check_output(idx, test): - # get results from listing file - fpth = os.path.join(test.workspace, f"{os.path.basename(test.name)}.lst") - budl = flopy.utils.Mf6ListBudget(fpth) - names = list(bud_lst) - d0 = budl.get_budget(names=names)[0] - dtype = d0.dtype - nbud = d0.shape[0] - - # get results from cbc file - cbc_bud = ["STO-SS", "STO-SY"] - d = np.recarray(nbud, dtype=dtype) - for key in bud_lst: - d[key] = 0.0 - fpth = os.path.join(test.workspace, f"{os.path.basename(test.name)}.cbc") - cobj = flopy.utils.CellBudgetFile(fpth, precision="double") - kk = cobj.get_kstpkper() - times = cobj.get_times() - for i, (k, t) in enumerate(zip(kk, times)): - for text in cbc_bud: - qin = 0.0 - qout = 0.0 - v = cobj.get_data(kstpkper=k, text=text)[0] - if isinstance(v, np.recarray): - vt = np.zeros(size3d, dtype=float) - for jdx, node in enumerate(v["node"]): - vt[node - 1] += v["q"][jdx] - v = vt.reshape(shape3d) - for kk in range(v.shape[0]): - for ii in range(v.shape[1]): - for jj in range(v.shape[2]): - vv = v[kk, ii, jj] - if vv < 0.0: - qout -= vv - else: - qin += vv - d["totim"][i] = t - d["time_step"][i] = k[0] - d["stress_period"] = k[1] - key = f"{text}_IN" - d[key][i] = qin - key = f"{text}_OUT" - d[key][i] = qout + from test_gwf_sto01 import check_output as check - diff = np.zeros((nbud, len(bud_lst)), dtype=float) - for i, key in enumerate(bud_lst): - diff[:, i] = d0[key] - d[key] - diffmax = np.abs(diff).max() - msg = f"maximum absolute total-budget difference ({diffmax}) " - - # write summary - fpth = os.path.join( - test.workspace, f"{os.path.basename(test.name)}.bud.cmp.out" - ) - with open(fpth, "w") as f: - for i in range(diff.shape[0]): - if i == 0: - line = f"{'TIME':>10s}" - for key in bud_lst: - line += f"{key + '_LST':>25s}" - line += f"{key + '_CBC':>25s}" - line += f"{key + '_DIF':>25s}" - f.write(line + "\n") - line = f"{d['totim'][i]:10g}" - for ii, key in enumerate(bud_lst): - line += f"{d0[key][i]:25g}" - line += f"{d[key][i]:25g}" - line += f"{diff[i, ii]:25g}" - f.write(line + "\n") - - if diffmax > budtol: - test.success = False - msg += f"exceeds {dtol}" - assert diffmax < dtol, msg - else: - test.success = True - print(" " + msg) + check(idx, test) # Check NetCDF output nc_fname = f"{os.path.basename(test.name)}.nc" @@ -397,6 +78,13 @@ def check_output(idx, test): ) hds = flopy.utils.HeadFile(hds_fpth, precision="double") + gwf = test.sims[0].gwf[0] + dis = getattr(gwf, "dis") + tdis = getattr(test.sims[0], "tdis") + nper = getattr(tdis, "nper").data + nstp = [1] + [6 for _ in range(nper - 1)] + nlay = getattr(dis, "nlay").data + # Compare NetCDF head arrays with binary headfile timestep = 0 for i in range(nper): diff --git a/autotest/test_netcdf_gwt_dsp01.py b/autotest/test_netcdf_gwt_dsp01.py index 96da1fd8125..24b3ab6f97e 100644 --- a/autotest/test_netcdf_gwt_dsp01.py +++ b/autotest/test_netcdf_gwt_dsp01.py @@ -1,3 +1,7 @@ +""" +NetCDF export test version of test_gwt_dsp01. +""" + import os import flopy @@ -12,214 +16,22 @@ except ImportError: pytest.skip("xuarray and xugrid not found", allow_module_level=True) -cases = ["dsp01a_ncf", "dsp01b_ncf"] -xt3d = [False, True] +cases = ["dsp01a", "dsp01b"] def build_models(idx, test): - nlay, nrow, ncol = 1, 1, 100 - nper = 1 - perlen = [5.0] - nstp = [200] - tsmult = [1.0] - steady = [True] - delr = 1.0 - delc = 1.0 - top = 1.0 - laytyp = 0 - ss = 0.0 - sy = 0.1 - botm = [0.0] - strt = 1.0 - hnoflo = 1e30 - hdry = -1e30 - hk = 1.0 - - c = {0: [[(0, 0, 0), 0.0000000], [(0, 0, 99), 0.0000000]]} - - nouter, ninner = 100, 300 - hclose, rclose, relax = 1e-6, 1e-6, 1.0 - - tdis_rc = [] - for i in range(nper): - tdis_rc.append((perlen[i], nstp[i], tsmult[i])) - - name = cases[idx] - - # build MODFLOW 6 files - ws = test.workspace - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws - ) - # create tdis package - tdis = flopy.mf6.ModflowTdis( - sim, - time_units="DAYS", - start_date_time="2041-01-01T00:00:00-05:00", - nper=nper, - perioddata=tdis_rc, - ) - - # create gwf model - gwfname = "gwf_" + name - gwf = flopy.mf6.MFModel( - sim, - model_type="gwf6", - modelname=gwfname, - model_nam_file=f"{gwfname}.nam", - ) - gwf.name_file.save_flows = True - - # create iterative model solution and register the gwf model with it - 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=f"{gwfname}.ims", - ) - sim.register_ims_package(imsgwf, [gwf.name]) + from test_gwt_dsp01 import build_models as build - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - idomain=np.ones((nlay, nrow, ncol), dtype=int), - filename=f"{gwfname}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{gwfname}.ic") - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, save_specific_discharge=True, icelltype=laytyp, k=hk, k33=hk - ) - - # chd files - chd = flopy.mf6.ModflowGwfchd( - gwf, - maxbound=len(c), - stress_period_data=c, - save_flows=False, - print_flows=True, - pname="CHD-1", - ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{gwfname}.cbc", - head_filerecord=f"{gwfname}.hds", - headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], - saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], - printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], - ) - - # create gwt model - gwtname = "gwt_" + name - gwt = flopy.mf6.MFModel( - sim, - model_type="gwt6", - modelname=gwtname, - model_nam_file=f"{gwtname}.nam", - ) - gwt.name_file.save_flows = True + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwt = sim.gwt[0] gwt.name_file.export_netcdf = "ugrid" - - # create iterative model solution and register the gwt model with it - 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=f"{gwtname}.ims", - ) - sim.register_ims_package(imsgwt, [gwt.name]) - - dis = flopy.mf6.ModflowGwtdis( - gwt, - export_array_netcdf=True, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - idomain=1, - filename=f"{gwtname}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwtic( - gwt, export_array_netcdf=True, strt=0.0, filename=f"{gwtname}.ic" - ) - - # advection - adv = flopy.mf6.ModflowGwtadv( - gwt, scheme="UPSTREAM", filename=f"{gwtname}.adv" - ) - - # dispersion - xt3d_off = not xt3d[idx] - dsp = flopy.mf6.ModflowGwtdsp( - gwt, - export_array_netcdf=True, - xt3d_off=xt3d_off, - diffc=100.0, - alh=0.0, - alv=0.0, - ath1=0.0, - atv=0.0, - filename=f"{gwtname}.dsp", - ) - - # constant concentration - fname = gwtname + ".cnc.obs.csv" - cnc_obs = {(fname): [("cnc000", "cnc", (0, 0, 0))]} - cnc_obs["digits"] = 15 - cncs = {0: [[(0, 0, 0), 1.0]]} - cnc = flopy.mf6.ModflowGwtcnc( - gwt, - maxbound=len(cncs), - stress_period_data=cncs, - print_flows=True, - save_flows=False, - observations=cnc_obs, - pname="CNC-1", - ) - - # mass storage and transfer - mst = flopy.mf6.ModflowGwtmst(gwt, porosity=0.1) - - # sources - ssm = flopy.mf6.ModflowGwtssm( - gwt, sources=[[]], print_flows=True, filename=f"{gwtname}.ssm" - ) + gwt.dis.export_array_netcdf = True + gwt.ic.export_array_netcdf = True + gwt.dsp.export_array_netcdf = True # output control + gwtname = "gwt_" + cases[idx] oc = flopy.mf6.ModflowGwtoc( gwt, budget_filerecord=f"{gwtname}.cbc", @@ -227,37 +39,17 @@ def build_models(idx, test): concentrationprintrecord=[ ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") ], - # saverecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], - # printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], - saverecord=[("CONCENTRATION", "ALL")], - printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "LAST")], + printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "LAST")], ) + return sim, dummy - # observations - obs_data0 = [("flow1", "flow-ja-face", (0, 0, 0), (0, 0, 1))] - obs_recarray = {f"{gwtname}.obs.csv": obs_data0} - obs = flopy.mf6.ModflowUtlobs( - gwt, - pname="gwt_obs", - filename=f"{gwtname}.obs", - digits=15, - print_input=True, - continuous=obs_recarray, - ) - # GWF GWT exchange - gwfgwt = flopy.mf6.ModflowGwfgwt( - sim, - exgtype="GWF6-GWT6", - exgmnamea=gwfname, - exgmnameb=gwtname, - filename=f"{name}.gwfgwt", - ) - - return sim, None +def check_output(idx, test): + from test_gwt_dsp01 import check_output as check + check(idx, test) -def check_output(idx, test): name = cases[idx] gwtname = "gwt_" + name @@ -270,133 +62,6 @@ def check_output(idx, test): except: assert False, f'could not load data from "{fpth}"' - # This is the answer to this problem. These concentrations are for - # time step 200. - cres = [ - [ - [ - 1.0, - 0.97472443, - 0.94947431, - 0.92427504, - 0.89915185, - 0.87412972, - 0.84923335, - 0.82448706, - 0.79991471, - 0.77553964, - 0.75138462, - 0.72747174, - 0.70382241, - 0.68045725, - 0.65739608, - 0.63465784, - 0.61226053, - 0.59022124, - 0.56855604, - 0.54727998, - 0.52640705, - 0.50595018, - 0.4859212, - 0.46633085, - 0.44718873, - 0.42850336, - 0.41028211, - 0.39253126, - 0.37525599, - 0.35846038, - 0.34214746, - 0.32631921, - 0.31097658, - 0.29611954, - 0.28174707, - 0.26785727, - 0.2544473, - 0.24151351, - 0.22905142, - 0.21705579, - 0.20552066, - 0.19443937, - 0.18380466, - 0.17360869, - 0.16384304, - 0.15449886, - 0.14556682, - 0.13703721, - 0.12889996, - 0.12114473, - 0.1137609, - 0.10673763, - 0.10006394, - 0.09372869, - 0.08772068, - 0.08202862, - 0.07664126, - 0.07154731, - 0.06673558, - 0.06219493, - 0.05791434, - 0.05388294, - 0.05009, - 0.04652499, - 0.04317758, - 0.04003765, - 0.03709534, - 0.03434103, - 0.03176537, - 0.0293593, - 0.02711402, - 0.02502107, - 0.02307224, - 0.02125967, - 0.01957579, - 0.01801336, - 0.01656542, - 0.01522538, - 0.01398691, - 0.01284404, - 0.01179109, - 0.01082267, - 0.00993375, - 0.00911954, - 0.0083756, - 0.00769775, - 0.00708212, - 0.00652511, - 0.00602341, - 0.00557398, - 0.00517407, - 0.00482116, - 0.00451303, - 0.0042477, - 0.00402344, - 0.00383879, - 0.00369253, - 0.00358368, - 0.00351152, - 0.00347556, - ] - ] - ] - cres = np.array(cres) - assert np.allclose( - cres, conc - ), "simulated concentrations do not match with known solution." - - # load the gwt observation file - fname = gwtname + ".obs.csv" - fname = os.path.join(test.workspace, fname) - gwtobs = np.genfromtxt(fname, names=True, delimiter=",", deletechars="") - - # load the cnc observation file - fname = gwtname + ".cnc.obs.csv" - fname = os.path.join(test.workspace, fname) - cncobs = np.genfromtxt(fname, names=True, delimiter=",", deletechars="") - - # ensure flow right face for first cell is equal to cnc flows - errmsg = f"observations not equal:\n{gwtobs}\n{cncobs}" - assert np.allclose(gwtobs["FLOW1"], -cncobs["CNC000"]), errmsg - # netcdf nc_fpth = os.path.join(test.workspace, f"{gwtname}.nc") ds = xu.open_dataset(nc_fpth) diff --git a/autotest/test_netcdf_gwt_prudic2004t2.py b/autotest/test_netcdf_gwt_prudic2004t2.py index b8db8def9ff..e6d3a11aecb 100644 --- a/autotest/test_netcdf_gwt_prudic2004t2.py +++ b/autotest/test_netcdf_gwt_prudic2004t2.py @@ -1,9 +1,5 @@ """ -Second problem described by Prudic et al (2004) -This problem involves transport through an aquifers, lakes and streams. -It requires the use of the Water Mover Package to send water from a stream, -into a lake, and then back into another stream. Solute is also transport -through the system. +NetCDF export test version of test_gwt_prudic2004t2. """ import os @@ -23,964 +19,31 @@ pytest.skip("xuarray and xugrid not found", allow_module_level=True) cases = ["prudic2004t2"] -data_path = project_root_path / "autotest" / "data" -model_path = data_path / "prudic2004test2" -fname = str(model_path / "lakibd.dat") -lakibd = np.loadtxt(fname, dtype=int) def build_models(idx, test): - ws = test.workspace - name = cases[idx] - gwfname = "gwf_" + name - gwtname = "gwt_" + name - sim = flopy.mf6.MFSimulation( - sim_name=name, - version="mf6", - exe_name="mf6", - sim_ws=ws, - continue_=False, - memory_print_option="ALL", - ) - - # number of time steps for period 2 are reduced from 12 * 25 to 25 in - # order to speed up this autotest - tdis_rc = [(1.0, 1, 1.0), (365.25 * 25, 25, 1.0)] - nper = len(tdis_rc) - tdis = flopy.mf6.ModflowTdis( - sim, - time_units="DAYS", - start_date_time="2041-01-01T00:00:00-05:00", - nper=nper, - perioddata=tdis_rc, - ) - - gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname) - - # ims - hclose = 0.001 - rclose = 0.1 - nouter = 1000 - ninner = 100 - relax = 0.99 - imsgwf = flopy.mf6.ModflowIms( - sim, - print_option="ALL", - 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=f"{gwfname}.ims", - ) - - nlay = 8 - nrow = 36 - ncol = 23 - delr = 405.665 - delc = 403.717 - top = 100.0 - fname = str(model_path / "bot1.dat") - bot0 = np.loadtxt(fname) - botm = [bot0] + [bot0 - (15.0 * k) for k in range(1, nlay)] - fname = str(model_path / "idomain1.dat") - idomain0 = np.loadtxt(fname, dtype=int) - idomain = nlay * [idomain0] - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - idomain=idomain, - ) - idomain = dis.idomain.array - - ic = flopy.mf6.ModflowGwfic(gwf, strt=50.0) - - npf = flopy.mf6.ModflowGwfnpf( - gwf, - xt3doptions=False, - save_flows=True, - save_specific_discharge=True, - icelltype=[1] + 7 * [0], - k=250.0, - k33=125.0, - ) - - sto_on = False - if sto_on: - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=True, - iconvert=[1] + 7 * [0], - ss=1.0e-5, - sy=0.3, - steady_state={0: True}, - transient={1: False}, - ) - - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{gwfname}.bud", - head_filerecord=f"{gwfname}.hds", - headprintrecord=[ - ("COLUMNS", ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - ) - - rch_on = True - if rch_on: - rch = flopy.mf6.ModflowGwfrcha( - gwf, recharge={0: 4.79e-3}, pname="RCH-1" - ) - - chdlist = [] - fname = str(model_path / "chd.dat") - for line in open(fname, "r").readlines(): - ll = line.strip().split() - if len(ll) == 4: - k, i, j, hd = ll - chdlist.append( - [ - ( - int(k) - 1, - int(i) - 1, - int(j) - 1, - ), - float(hd), - ] - ) - chd = flopy.mf6.ModflowGwfchd( - gwf, stress_period_data=chdlist, pname="CHD-1" - ) - - rivlist = [] - fname = str(model_path / "riv.dat") - for line in open(fname, "r").readlines(): - ll = line.strip().split() - if len(ll) == 7: - k, i, j, s, c, rb, bn = ll - rivlist.append( - [ - ( - int(k) - 1, - int(i) - 1, - int(j) - 1, - ), - float(s), - float(c), - float(rb), - bn, - ] - ) - rivra = flopy.mf6.ModflowGwfriv.stress_period_data.empty( - gwf, maxbound=len(rivlist), boundnames=True - )[0] - for i, t in enumerate(rivlist): - rivra[i] = tuple(t) - sfrpd = np.genfromtxt(model_path / "sfr-packagedata.dat", names=True) - sfrpackagedata = flopy.mf6.ModflowGwfsfr.packagedata.empty( - gwf, boundnames=True, maxbound=sfrpd.shape[0] - ) - for name in sfrpackagedata.dtype.names: - if name in rivra.dtype.names: - sfrpackagedata[name] = rivra[name] - for name in sfrpackagedata.dtype.names: - if name in sfrpd.dtype.names: - sfrpackagedata[name] = sfrpd[name] - sfrpackagedata["boundname"] = rivra["boundname"] - with open(model_path / "sfr-connectiondata.dat") as f: - lines = f.readlines() - sfrconnectiondata = [] - for line in lines: - t = line.split() - c = [] - for v in t: - i = int(v) - c.append(i) - sfrconnectiondata.append(c) - sfrperioddata = {0: [[0, "inflow", 86400], [18, "inflow", 8640.0]]} - - sfr_on = True - if sfr_on: - sfr = flopy.mf6.ModflowGwfsfr( - gwf, - print_stage=True, - print_flows=True, - stage_filerecord=gwfname + ".sfr.bin", - budget_filerecord=gwfname + ".sfr.bud", - mover=True, - pname="SFR-1", - length_conversion=3.28084, - time_conversion=86400.0, - boundnames=True, - nreaches=len(rivlist), - packagedata=sfrpackagedata, - connectiondata=sfrconnectiondata, - perioddata=sfrperioddata, - ) - - lakeconnectiondata = [] - nlakecon = [0, 0] - lak_leakance = 1.0 - for i in range(nrow): - for j in range(ncol): - if lakibd[i, j] == 0: - continue - else: - ilak = lakibd[i, j] - 1 - # back - if i > 0: - if lakibd[i - 1, j] == 0 and idomain[0, i - 1, j]: - h = [ - ilak, - nlakecon[ilak], - (0, i - 1, j), - "horizontal", - lak_leakance, - 0.0, - 0.0, - delc / 2.0, - delr, - ] - nlakecon[ilak] += 1 - lakeconnectiondata.append(h) - # left - if j > 0: - if lakibd[i, j - 1] and idomain[0, i, j - 1] == 0: - h = [ - ilak, - nlakecon[ilak], - (0, i, j - 1), - "horizontal", - lak_leakance, - 0.0, - 0.0, - delr / 2.0, - delc, - ] - nlakecon[ilak] += 1 - lakeconnectiondata.append(h) - # right - if j < ncol - 1: - if lakibd[i, j + 1] == 0 and idomain[0, i, j + 1]: - h = [ - ilak, - nlakecon[ilak], - (0, i, j + 1), - "horizontal", - lak_leakance, - 0.0, - 0.0, - delr / 2.0, - delc, - ] - nlakecon[ilak] += 1 - lakeconnectiondata.append(h) - # front - if i < nrow - 1: - if lakibd[i + 1, j] == 0 and idomain[0, i + 1, j]: - h = [ - ilak, - nlakecon[ilak], - (0, i + 1, j), - "horizontal", - lak_leakance, - 0.0, - 0.0, - delc / 2.0, - delr, - ] - nlakecon[ilak] += 1 - lakeconnectiondata.append(h) - # vertical - v = [ - ilak, - nlakecon[ilak], - (1, i, j), - "vertical", - lak_leakance, - 0.0, - 0.0, - 0.0, - 0.0, - ] - nlakecon[ilak] += 1 - lakeconnectiondata.append(v) - - i, j = np.where(lakibd > 0) - idomain[0, i, j] = 0 - gwf.dis.idomain.set_data(idomain[0], layer=0, multiplier=[1]) - - lakpackagedata = [ - [0, 44.0, nlakecon[0], "lake1"], - [1, 35.2, nlakecon[1], "lake2"], - ] - # - outlets = [[0, 0, -1, "MANNING", 44.5, 5.000000, 0.03, 0.2187500e-02]] - - lake_on = True - if lake_on: - lak = flopy.mf6.ModflowGwflak( - gwf, - time_conversion=86400.000, - print_stage=True, - print_flows=True, - stage_filerecord=gwfname + ".lak.bin", - budget_filerecord=gwfname + ".lak.bud", - mover=True, - pname="LAK-1", - boundnames=True, - nlakes=2, - noutlets=len(outlets), - outlets=outlets, - packagedata=lakpackagedata, - connectiondata=lakeconnectiondata, - ) - - mover_on = True - if mover_on: - maxmvr, maxpackages = 2, 2 - mvrpack = [["SFR-1"], ["LAK-1"]] - mvrperioddata = [ - ["SFR-1", 5, "LAK-1", 0, "FACTOR", 1.0], - ["LAK-1", 0, "SFR-1", 6, "FACTOR", 1.0], - ] - mvr = flopy.mf6.ModflowGwfmvr( - gwf, - maxmvr=maxmvr, - print_flows=True, - maxpackages=maxpackages, - packages=mvrpack, - perioddata=mvrperioddata, - ) - - transport = True - if transport: - gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname) - gwt.name_file.export_netcdf = "ugrid" - - # ims - hclose = 0.001 - rclose = 0.001 - nouter = 50 - ninner = 20 - relax = 0.97 - imsgwt = flopy.mf6.ModflowIms( - sim, - print_option="ALL", - outer_dvclose=hclose, - outer_maximum=nouter, - inner_maximum=ninner, - under_relaxation="DBD", - under_relaxation_theta=0.7, - inner_dvclose=hclose, - rcloserecord=rclose, - linear_acceleration="BICGSTAB", - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=relax, - filename=f"{gwtname}.ims", - ) - sim.register_ims_package(imsgwt, [gwt.name]) - - dis = flopy.mf6.ModflowGwtdis( - gwt, - export_array_netcdf=True, - nlay=nlay, - nrow=nrow, - ncol=ncol, - delr=delr, - delc=delc, - top=top, - botm=botm, - idomain=idomain, - ) - ic = flopy.mf6.ModflowGwtic(gwt, export_array_netcdf=True, strt=0.0) - mst = flopy.mf6.ModflowGwtmst(gwt, porosity=0.3) - adv = flopy.mf6.ModflowGwtadv(gwt, scheme="TVD") - dsp = flopy.mf6.ModflowGwtdsp( - gwt, export_array_netcdf=True, alh=20.0, ath1=2, atv=0.2 - ) - sourcerecarray = [()] - ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) - cnclist = [ - [(0, 0, 11), 500.0], - [(0, 0, 12), 500.0], - [(0, 0, 13), 500.0], - [(0, 0, 14), 500.0], - [(1, 0, 11), 500.0], - [(1, 0, 12), 500.0], - [(1, 0, 13), 500.0], - [(1, 0, 14), 500.0], - ] - cnc = flopy.mf6.ModflowGwtcnc( - gwt, - maxbound=len(cnclist), - stress_period_data=cnclist, - save_flows=False, - pname="CNC-1", - ) - - lktpackagedata = [ - (0, 0.0, 99.0, 999.0, "mylake1"), - (1, 0.0, 99.0, 999.0, "mylake2"), - ] - lktperioddata = [ - (0, "STATUS", "ACTIVE"), - (1, "STATUS", "ACTIVE"), - ] - - lkt_obs = {} - for obstype in [ - "CONCENTRATION", - "STORAGE", - "CONSTANT", - "FROM-MVR", - "TO-MVR", - "RAINFALL", - "EVAPORATION", - "RUNOFF", - "EXT-INFLOW", - "WITHDRAWAL", - "EXT-OUTFLOW", - ]: - fname = f"{gwtname}.lkt.obs.{obstype.lower()}.csv" - obs1 = [] - ncv = 2 - if obstype == "TO-MVR": - ncv = 1 - obs1 = [(f"lkt{i + 1}", obstype, i + 1) for i in range(ncv)] - obs2 = [ - (f"blkt{i + 1}", obstype, f"mylake{i + 1}") for i in range(ncv) - ] - lkt_obs[fname] = obs1 + obs2 - - # add LKT specific obs - obstype = "LKT" - ncv = 2 - nconn = [67, 32] - fname = f"{gwtname}.lkt.obs.{obstype.lower()}.csv" - obs1 = [ - (f"lkt{1}-{iconn + 1}", obstype, 1, iconn + 1) - for iconn in range(nconn[0]) # lake 1 - ] + [ - (f"lkt{2}-{iconn + 1}", obstype, 2, iconn + 1) - for iconn in range(nconn[1]) # lake 2 - ] - obs2 = [ - (f"blkt{i + 1}", obstype, f"mylake{i + 1}") for i in range(ncv) - ] - lkt_obs[fname] = obs1 + obs2 - - lkt_obs["digits"] = 15 - lkt_obs["print_input"] = True - lkt_obs["filename"] = gwtname + ".lkt.obs" - - lkt_on = True - if lkt_on: - lkt = flopy.mf6.modflow.ModflowGwtlkt( - gwt, - boundnames=True, - save_flows=True, - print_input=True, - print_flows=True, - print_concentration=True, - concentration_filerecord=gwtname + ".lkt.bin", - budget_filerecord=gwtname + ".lkt.bud", - packagedata=lktpackagedata, - lakeperioddata=lktperioddata, - observations=lkt_obs, - pname="LAK-1", - auxiliary=["aux1", "aux2"], - ) - - sftpackagedata = [] - for irno in range(sfrpd.shape[0]): - t = (irno, 0.0, 99.0, 999.0, f"myreach{irno + 1}") - sftpackagedata.append(t) - - sftperioddata = [(0, "STATUS", "ACTIVE"), (0, "CONCENTRATION", 0.0)] - - sft_obs = {} - for obstype in [ - "CONCENTRATION", - "STORAGE", - "CONSTANT", - "FROM-MVR", - "TO-MVR", - "SFT", - "RAINFALL", - "EVAPORATION", - "RUNOFF", - "EXT-INFLOW", - "EXT-OUTFLOW", - ]: - fname = f"{gwtname}.sft.obs.{obstype.lower()}.csv" - obs1 = [ - (f"sft{i + 1}", obstype, i + 1) for i in range(sfrpd.shape[0]) - ] - obs2 = [ - (f"bsft{i + 1}", obstype, f"myreach{i + 1}") - for i in range(sfrpd.shape[0]) - ] - sft_obs[fname] = obs1 + obs2 - - obstype = "FLOW-JA-FACE" - fname = f"{gwtname}.sft.obs.{obstype.lower()}.csv" - obs1 = [] - for id1, reach_connections in enumerate(sfrconnectiondata): - for id2 in reach_connections: - id2 = abs(id2) - if id1 != id2: - obs1.append( - (f"sft{id1 + 1}x{id2 + 1}", obstype, id1 + 1, id2 + 1) - ) - obs2 = [ - (f"bsft{i + 1}", obstype, f"myreach{i + 1}") - for i in range(sfrpd.shape[0]) - ] - sft_obs[fname] = obs1 + obs2 - - # append additional obs attributes to obs dictionary - sft_obs["digits"] = 15 - sft_obs["print_input"] = True - sft_obs["filename"] = gwtname + ".sft.obs" - - sft_on = True - if sft_on: - sft = flopy.mf6.modflow.ModflowGwtsft( - gwt, - boundnames=True, - save_flows=True, - print_input=True, - print_flows=True, - print_concentration=True, - concentration_filerecord=gwtname + ".sft.bin", - budget_filerecord=gwtname + ".sft.bud", - packagedata=sftpackagedata, - reachperioddata=sftperioddata, - observations=sft_obs, - pname="SFR-1", - auxiliary=["aux1", "aux2"], - ) - - # mover transport package - mvt = flopy.mf6.modflow.ModflowGwtmvt(gwt) - - oc = flopy.mf6.ModflowGwtoc( - gwt, - budget_filerecord=f"{gwtname}.cbc", - concentration_filerecord=f"{gwtname}.ucn", - concentrationprintrecord=[ - ("COLUMNS", ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], - printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], - ) - - # GWF GWT exchange - gwfgwt = flopy.mf6.ModflowGwfgwt( - sim, exgtype="GWF6-GWT6", exgmnamea=gwfname, exgmnameb=gwtname - ) - - return sim, None - + from test_gwt_prudic2004t2 import build_models as build -def make_concentration_vs_time(sim): - print("making plot of concentration versus time...") - name = sim.name - ws = sim.workspace - sim = flopy.mf6.MFSimulation.load(sim_ws=ws) - gwfname = "gwf_" + name - gwtname = "gwt_" + name - gwf = sim.get_model(gwfname) - gwt = sim.get_model(gwtname) - - bobj = gwt.lkt.output.concentration() - lkaconc = bobj.get_alldata()[:, 0, 0, :] - times = bobj.times - bobj.file.close() - - bobj = gwt.sft.output.concentration() - sfaconc = bobj.get_alldata()[:, 0, 0, :] - times = bobj.times - bobj.file.close() - - import matplotlib.pyplot as plt - - plt.figure(figsize=(8, 5)) - times = np.array(times) / 365.0 - plt.plot(times, lkaconc[:, 0], "b-", label="Lake 1") - plt.plot(times, sfaconc[:, 30], "r-", label="Stream segment 3") - plt.plot(times, sfaconc[:, 37], "g-", label="Stream segment 4") - plt.legend() - plt.ylim(0, 50) - plt.xlim(0, 25) - plt.xlabel("TIME, IN YEARS") - plt.ylabel("SIMULATED BORON CONCENTRATION,\nIN MICROGRAMS PER LITER") - plt.draw() - fname = os.path.join(ws, "fig-concentration_vs_time.png") - print(f"Creating {fname}") - plt.savefig(fname) - - return - - -def make_concentration_map(sim): - print("making concentration map...") - - import matplotlib.pyplot as plt - - levels = [ - 1, - 10, - 25, - 50, - 100, - 150, - 200, - 250, - 300, - 350, - 400, - 450, - 500, - ] - - name = sim.name - ws = sim.workspace - simfp = flopy.mf6.MFSimulation.load(sim_ws=ws) - gwfname = "gwf_" + name - gwtname = "gwt_" + name - gwf = simfp.get_model(gwfname) - gwt = simfp.get_model(gwtname) - conc = gwt.output.concentration().get_data() - lakconc = gwt.lak.output.concentration().get_data().flatten() - - il, jl = np.where(lakibd > 0) - for i, j in zip(il, jl): - ilak = lakibd[i, j] - 1 - lake_conc = lakconc[ilak] - conc[0, i, j] = lake_conc - - fig, axs = plt.subplots(2, 2, figsize=(5, 7), dpi=300, tight_layout=True) - - # plot layers 1, 3, 5, and 8 - for iplot, ilay in enumerate([0, 2, 4, 7]): - ax = axs.flatten()[iplot] - pmv = flopy.plot.PlotMapView(model=gwt, ax=ax, layer=ilay) - # pmv.plot_grid() - pmv.plot_array(lakibd, masked_values=[0], alpha=0.2) - pmv.plot_inactive(color_noflow="gray", alpha=0.25) - # pmv.plot_bc(name="CHD-1", color="blue") - cs = pmv.contour_array(conc, levels=levels, masked_values=[1.0e30]) - ax.clabel(cs, cs.levels[::1], fmt="%1.0f", colors="b") - ax.set_title(f"Model layer {ilay + 1}") - - fname = os.path.join(ws, "fig-concentration.png") - print(f"Creating {fname}") - plt.savefig(fname) - - return - - -def check_obs(sim): - print("checking obs...") - name = sim.name - ws = sim.workspace - sim = flopy.mf6.MFSimulation.load(sim_ws=ws) - gwfname = "gwf_" + name - gwtname = "gwt_" + name - gwf = sim.get_model(gwfname) - gwt = sim.get_model(gwtname) - - # ensure SFT obs are the same whether specified by - # boundname or by reach - csvfiles = gwt.sft.obs.output.obs_names - for csvfile in csvfiles: - if ".flow-ja-face.csv" in csvfile: - continue - print(f"Checking csv file: {csvfile}") - conc_ra = gwt.sft.obs.output.obs(f=csvfile).data - # save a couple entries for comparison with lake - if ".to-mvr." in csvfile: - sft6tomvr = conc_ra[f"BSFT6"] - if ".from-mvr." in csvfile: - sft7tomvr = conc_ra[f"BSFT7"] - success = True - for ireach in range(38): - # print(f" Checking reach {ireach + 1}") - is_same = np.allclose( - conc_ra[f"SFT{ireach + 1}"], conc_ra[f"BSFT{ireach + 1}"] - ) - if not is_same: - success = False - for t, x, y in zip( - conc_ra["totim"], - conc_ra[f"SFT{ireach + 1}"], - conc_ra[f"BSFT{ireach + 1}"], - ): - print(t, x, y) - - # process the sft values and make sure the individual connection rates add up to the boundname rate - csvfile = "gwt_prudic2004t2.sft.obs.flow-ja-face.csv" - print(f"Checking csv file: {csvfile}") - conc_ra = gwt.sft.obs.output.obs(f=csvfile).data - ntimes = conc_ra.shape[0] - for ireach in range(38): - connection_sum = np.zeros(ntimes) - for column_name in conc_ra.dtype.names: - if f"SFT{ireach + 1}X" in column_name: - connection_sum += conc_ra[column_name] - is_same = np.allclose(connection_sum, conc_ra[f"BSFT{ireach + 1}"]) - if not is_same: - success = False - diff = connection_sum - conc_ra[f"BSFT{ireach + 1}"] - print( - f"Problem with SFT {ireach + 1}; mindiff {diff.min()} and maxdiff {diff.max()}" - ) - # for itime, (cs, bsft) in enumerate(zip(connection_sum, conc_ra[f"BSFT{ireach + 1}"])): - # print(itime, cs, bsft) - - assert success, "One or more SFT obs checks did not pass" - - # ensure LKT obs are the same whether specified by - # boundname or by reach - csvfiles = gwt.lkt.obs.output.obs_names - for csvfile in csvfiles: - if ".lkt.csv" in csvfile: - continue - print(f"Checking csv file: {csvfile}") - conc_ra = gwt.lkt.obs.output.obs(f=csvfile).data - if ".from-mvr." in csvfile: - lkt1frommvr = conc_ra[f"BLKT1"] - if ".to-mvr." in csvfile: - lkt1tomvr = conc_ra[f"BLKT1"] - success = True - if ".to-mvr." in csvfile: - numvalues = 1 # outlet - else: - numvalues = 2 # lakes - for ilake in range(numvalues): - # print(f" Checking lake {ilake + 1}") - is_same = np.allclose( - conc_ra[f"LKT{ilake + 1}"], conc_ra[f"BLKT{ilake + 1}"] - ) - if not is_same: - success = False - for t, x, y in zip( - conc_ra["totim"], - conc_ra[f"LKT{ilake + 1}"], - conc_ra[f"BLKT{ilake + 1}"], - ): - print(t, x, y) - - # process the lkt values and make sure the individual connection rates add up to the boundname rate - csvfile = "gwt_prudic2004t2.lkt.obs.lkt.csv" - print(f"Checking csv file: {csvfile}") - conc_ra = gwt.lkt.obs.output.obs(f=csvfile).data - ntimes = conc_ra.shape[0] - for ilake in [0, 1]: - connection_sum = np.zeros(ntimes) - for column_name in conc_ra.dtype.names: - if f"LKT{ilake + 1}" in column_name and column_name.startswith( - "LKT" - ): - connection_sum += conc_ra[column_name] - is_same = np.allclose(connection_sum, conc_ra[f"BLKT{ilake + 1}"]) - if not is_same: - success = False - print(f"Problem with Lake {ilake + 1}") - for itime, (cs, blkt) in enumerate( - zip(connection_sum, conc_ra[f"BLKT1"]) - ): - print(itime, cs, blkt) - - assert success, "One or more LKT obs checks did not pass" - - # check that SFT6 to-mvr is equal to LKT1 from-mvr - success = True - is_same = np.allclose(-sft6tomvr, lkt1frommvr, atol=0.1) - if not is_same: - success = False - print(f"Problem with sft6tomvr comparison to lkt1frommvr") - for itime, (a, b) in enumerate(zip(-sft6tomvr, lkt1frommvr)): - print(itime, a, b) - - is_same = np.allclose(-lkt1tomvr, sft7tomvr) - if not is_same: - success = False - print(f"Problem with lkt1tomvr comparison to sft7tomvr") - for itime, (a, b) in enumerate(zip(-lkt1tomvr, sft7tomvr)): - print(itime, a, b) - - assert success, "One or more SFT-LKT obs checks did not pass" + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwt = sim.gwt[0] + gwt.name_file.export_netcdf = "ugrid" + gwt.dis.export_array_netcdf = True + gwt.ic.export_array_netcdf = True + gwt.dsp.export_array_netcdf = True + return sim, dummy def check_output(idx, test): - makeplot = False - for arg in sys.argv: - if arg.lower() == "--makeplot": - makeplot = True + from test_gwt_prudic2004t2 import check_output as check - if makeplot: - make_concentration_vs_time(test) - make_concentration_map(test) + check(idx, test) - # ensure concentrations were saved + # netcdf ws = test.workspace name = test.name gwtname = "gwt_" + name - - check_obs(test) - - fname = gwtname + ".lkt.bin" - fname = os.path.join(ws, fname) - bobj = flopy.utils.HeadFile( - fname, precision="double", text="concentration" - ) - lkaconc = bobj.get_alldata()[:, 0, 0, :] - times = bobj.times - bobj.file.close() - - fname = gwtname + ".sft.bin" - fname = os.path.join(ws, fname) - bobj = flopy.utils.HeadFile( - fname, precision="double", text="concentration" - ) - sfaconc = bobj.get_alldata()[:, 0, 0, :] - times = bobj.times - bobj.file.close() - - # set atol - atol = 0.02 - - # check simulated concentration in lak 1 and 2 sfr reaches - res_lak1 = lkaconc[:, 0] - ans_lak1 = [ - -1.73249951e-19, - 5.97398983e-02, - 4.18358112e-01, - 1.48857598e00, - 3.63202585e00, - 6.92925430e00, - 1.11162776e01, - 1.57328143e01, - 2.03088745e01, - 2.45013060e01, - 2.81200704e01, - 3.11132152e01, - 3.34833369e01, - 3.53028319e01, - 3.66693021e01, - 3.76781530e01, - 3.84188513e01, - 3.89615387e01, - 3.93577458e01, - 3.96464993e01, - 3.98598113e01, - 4.00184878e01, - 4.01377654e01, - 4.02288674e01, - 4.02998291e01, - 4.03563314e01, - ] - ans_lak1 = np.array(ans_lak1) - d = res_lak1 - ans_lak1 - msg = f"{res_lak1} {ans_lak1} {d}" - assert np.allclose(res_lak1, ans_lak1, atol=atol), msg - - res_sfr3 = sfaconc[:, 30] - ans_sfr3 = [ - -7.67944651e-23, - 5.11358249e-03, - 3.76169957e-02, - 1.42055634e-01, - 3.72438193e-01, - 7.74112522e-01, - 1.37336373e00, - 2.18151231e00, - 3.19993561e00, - 4.42853144e00, - 5.85660993e00, - 7.46619448e00, - 9.22646330e00, - 1.11069607e01, - 1.30764504e01, - 1.50977917e01, - 1.71329980e01, - 1.91636634e01, - 2.11530199e01, - 2.30688490e01, - 2.48821059e01, - 2.65691424e01, - 2.81080543e01, - 2.94838325e01, - 3.06909748e01, - 3.17352915e01, - ] - ans_sfr3 = np.array(ans_sfr3) - d = res_sfr3 - ans_sfr3 - msg = f"{res_sfr3} {ans_sfr3} {d}" - assert np.allclose(res_sfr3, ans_sfr3, atol=atol), msg - - res_sfr4 = sfaconc[:, 37] - ans_sfr4 = [ - -2.00171747e-20, - 3.55076535e-02, - 2.49465789e-01, - 8.91299656e-01, - 2.18622372e00, - 4.19920114e00, - 6.79501651e00, - 9.72255743e00, - 1.27208739e01, - 1.55989390e01, - 1.82462345e01, - 2.06258607e01, - 2.27255881e01, - 2.45721928e01, - 2.62061367e01, - 2.76640442e01, - 2.89788596e01, - 3.01814571e01, - 3.12842113e01, - 3.22945541e01, - 3.32174210e01, - 3.40539043e01, - 3.48027700e01, - 3.54636082e01, - 3.60384505e01, - 3.65330352e01, - ] - ans_sfr4 = np.array(ans_sfr4) - d = res_sfr4 - ans_sfr4 - msg = f"{res_sfr4} {ans_sfr4} {d}" - assert np.allclose(res_sfr4, ans_sfr4, atol=atol), msg - - # used to make results for the gwtgwt version of this problem - # fname = os.path.join(ws, f"result_conc_lak1.txt") - # np.savetxt(fname, res_lak1) - # fname = os.path.join(ws, f"result_conc_sfr3.txt") - # np.savetxt(fname, res_sfr3) - # fname = os.path.join(ws, f"result_conc_sfr4.txt") - # np.savetxt(fname, res_sfr4) - - # netcdf nc_fpth = os.path.join(ws, f"{gwtname}.nc") - print(nc_fpth) ds = xu.open_dataset(nc_fpth) xds = ds.ugrid.to_dataset() @@ -1004,8 +67,9 @@ def check_output(idx, test): for l in range(nlay): assert np.allclose( np.array(rec[l]).flatten(), - #xds[f"concentration_l{l+1}"][timestep, :].data, - xds[f"concentration_l{l+1}"][timestep, :].fillna(1.00000000e+30).data, + xds[f"concentration_l{l+1}"][timestep, :] + .fillna(1.00000000e30) + .data, ), f"NetCDF-concentration comparison failure in timestep {timestep+1}" timestep += 1 @@ -1049,6 +113,7 @@ def check_output(idx, test): @pytest.mark.slow +@pytest.mark.netcdf @pytest.mark.parametrize("idx, name", enumerate(cases)) def test_mf6model(idx, name, function_tmpdir, targets): test = TestFramework(