diff --git a/.docs/Notebooks/netcdf01_tutorial.py b/.docs/Notebooks/netcdf01_tutorial.py new file mode 100644 index 000000000..b70c22af3 --- /dev/null +++ b/.docs/Notebooks/netcdf01_tutorial.py @@ -0,0 +1,344 @@ +import sys +from pathlib import Path +from pprint import pformat +from tempfile import TemporaryDirectory + +import numpy as np +import xarray as xr + +import flopy + +print(sys.version) +print(f"flopy version: {flopy.__version__}") + +DNODATA = 3.0e30 + + +def create_sim(ws): + name = "uzf01" + perlen = [500.0] + nper = len(perlen) + nstp = [10] + tsmult = nper * [1.0] + crs = "EPSG:26916" + nlay, nrow, ncol = 100, 1, 1 + delr = 1.0 + delc = 1.0 + delv = 1.0 + top = 100.0 + botm = [top - (k + 1) * delv for k in range(nlay)] + strt = 0.5 + hk = 1.0 + laytyp = 1 + ss = 0.0 + sy = 0.1 + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + # build MODFLOW 6 files + 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", nper=nper, perioddata=tdis_rc) + + # create iterative model solution and register the gwf model with it + nouter, ninner = 100, 10 + hclose, rclose, relax = 1.5e-6, 1e-6, 0.97 + imsgwf = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.7, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + + # create gwf model + newtonoptions = "NEWTON UNDER_RELAXATION" + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name, + newtonoptions=newtonoptions, + save_flows=True, + ) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + crs=crs, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=np.ones((nlay, nrow, ncol), dtype=int), + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=False, icelltype=laytyp, k=hk) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=laytyp, + ss=ss, + sy=sy, + steady_state={0: False}, + transient={0: True}, + ) + + # ghbg + ghb_obs = {f"{name}.ghb.obs.csv": [("100_1_1", "GHB", (99, 0, 0))]} + bhead = np.full(nlay * nrow * ncol, DNODATA, dtype=float) + cond = np.full(nlay * nrow * ncol, DNODATA, dtype=float) + bhead[nlay - 1] = 1.5 + cond[nlay - 1] = 1.0 + ghb = flopy.mf6.ModflowGwfghbg( + gwf, + print_input=True, + print_flows=True, + bhead=bhead, + cond=cond, + save_flows=False, + ) + + ghb.obs.initialize( + filename=f"{name}.ghb.obs", + digits=20, + print_input=True, + continuous=ghb_obs, + ) + + # note: for specifying lake number, use fortran indexing! + uzf_obs = { + f"{name}.uzf.obs.csv": [ + ("wc 02", "water-content", 2, 0.5), + ("wc 50", "water-content", 50, 0.5), + ("wcbn 02", "water-content", "uzf 002", 0.5), + ("wcbn 50", "water-content", "UZF 050", 0.5), + ("rch 02", "uzf-gwrch", "uzf 002"), + ("rch 50", "uzf-gwrch", "uzf 050"), + ] + } + + sd = 0.1 + vks = hk + thtr = 0.05 + thti = thtr + thts = sy + eps = 4 + uzf_pkdat = [[0, (0, 0, 0), 1, 1, sd, vks, thtr, thts, thti, eps, "uzf 001"]] + [ + [k, (k, 0, 0), 0, k + 1, sd, vks, thtr, thts, thti, eps, f"uzf {k + 1:03d}"] + for k in range(1, nlay - 1) + ] + uzf_pkdat[-1][3] = -1 + infiltration = 2.01 + uzf_spd = {0: [[0, infiltration, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]} + uzf = flopy.mf6.ModflowGwfuzf( + gwf, + print_input=True, + print_flows=True, + save_flows=True, + boundnames=True, + ntrailwaves=15, + nwavesets=40, + nuzfcells=len(uzf_pkdat), + packagedata=uzf_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{name}.uzf.bud", + budgetcsv_filerecord=f"{name}.uzf.bud.csv", + observations=uzf_obs, + filename=f"{name}.uzf", + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.bud", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")], + ) + + obs_lst = [] + obs_lst.append(["obs1", "head", (0, 0, 0)]) + obs_lst.append(["obs2", "head", (1, 0, 0)]) + obs_dict = {f"{name}.obs.csv": obs_lst} + obs = flopy.mf6.ModflowUtlobs(gwf, pname="head_obs", digits=20, continuous=obs_dict) + + return sim + + +def add_netcdf_vars(dataset, nc_info, dimmap): + def _data_shape(shape): + dims_l = [] + for d in shape: + dims_l.append(dimmap[d]) + + return dims_l + + for v in nc_info: + varname = nc_info[v]["varname"] + data = np.full( + _data_shape(nc_info[v]["nc_shape"]), + nc_info[v]["attrs"]["_FillValue"], + dtype=nc_info[v]["nc_type"], + ) + var_d = {varname: (nc_info[v]["nc_shape"], data)} + dataset = dataset.assign(var_d) + for a in nc_info[v]["attrs"]: + dataset[varname].attrs[a] = nc_info[v]["attrs"][a] + + return dataset + + +temp_dir = TemporaryDirectory() +workspace = Path(temp_dir.name) + +# run the non-netcdf simulation +sim = create_sim(ws=workspace) +sim.write_simulation() +success, buff = sim.run_simulation(silent=True, report=True) +assert success, pformat(buff) + +# create directory for netcdf sim +sim.set_sim_path(workspace / "netcdf") +gwf = sim.get_model("uzf01") +gwf.name_file.nc_filerecord = "uzf01.structured.nc" +sim.write_simulation() + +# create the netcdf dataset +ds = xr.Dataset() + +# get model netcdf info +nc_info = gwf.netcdf_info() + +# update dataset with required attributes +for a in nc_info["attrs"]: + ds.attrs[a] = nc_info["attrs"][a] + +# set dimensional info +dis = gwf.modelgrid +xoff = dis.xoffset +yoff = dis.yoffset +x = xoff + dis.xycenters[0] +y = yoff + dis.xycenters[1] +z = [float(x) for x in range(1, dis.nlay + 1)] +nstp = sum(gwf.modeltime.nstp) +time = gwf.modeltime.tslen +nlay = dis.nlay +nrow = dis.nrow +ncol = dis.ncol +dimmap = {"time": nstp, "z": nlay, "y": nrow, "x": ncol} + +# create coordinate vars +var_d = {"time": (["time"], time), "z": (["z"], z), "y": (["y"], y), "x": (["x"], x)} +ds = ds.assign(var_d) + +# dis +dis = gwf.get_package("dis") +nc_info = dis.netcdf_info() + +# create dis dataset variables +ds = add_netcdf_vars(ds, nc_info, dimmap) + +# update data +ds["dis_delr"].values = dis.delr.get_data() +ds["dis_delc"].values = dis.delc.get_data() +ds["dis_top"].values = dis.top.get_data() +ds["dis_botm"].values = dis.botm.get_data() +ds["dis_idomain"].values = dis.idomain.get_data() + +# update dis to read from netcdf +with open(workspace / "netcdf" / "uzf01.dis", "w") as f: + f.write("BEGIN options\n") + f.write(" crs EPSG:26916\n") + f.write("END options\n\n") + f.write("BEGIN dimensions\n") + f.write(" NLAY 100\n") + f.write(" NROW 1\n") + f.write(" NCOL 1\n") + f.write("END dimensions\n\n") + f.write("BEGIN griddata\n") + f.write(" delr NETCDF\n") + f.write(" delc NETCDF\n") + f.write(" top NETCDF\n") + f.write(" botm NETCDF\n") + f.write(" idomain NETCDF\n") + f.write("END griddata\n") + +# npf +npf = gwf.get_package("npf") +nc_info = npf.netcdf_info() + +# create npf dataset variables +ds = add_netcdf_vars(ds, nc_info, dimmap) + +# update data +ds["npf_icelltype"].values = npf.icelltype.get_data() +ds["npf_k"].values = npf.k.get_data() + +# update npf to read from netcdf +with open(workspace / "netcdf" / "uzf01.npf", "w") as f: + f.write("BEGIN options\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" icelltype NETCDF\n") + f.write(" k NETCDF\n") + f.write("END griddata\n") + +# get ghbg package netcdf info +ghbg = gwf.get_package("ghbg_0") +nc_info = ghbg.netcdf_info() + +# create ghbg dataset variables +ds = add_netcdf_vars(ds, nc_info, dimmap) + +# update bhead netcdf array from flopy perioddata +for p in ghbg.bhead.get_data(): + istp = sum(gwf.modeltime.nstp[0:p]) + ds["ghbg_0_bhead"].values[istp] = ghbg.bhead.get_data()[p] + +# update cond netcdf array from flopy perioddata +for p in ghbg.cond.get_data(): + istp = sum(gwf.modeltime.nstp[0:p]) + ds["ghbg_0_cond"].values[istp] = ghbg.cond.get_data()[p] + +# write the netcdf +ds.to_netcdf( + workspace / "netcdf/uzf01.structured.nc", format="NETCDF4", engine="netcdf4" +) + +# update ghbg to read from netcdf +with open(workspace / "netcdf/uzf01.ghbg", "w") as f: + f.write("BEGIN options\n") + f.write(" READARRAYGRID\n") + f.write(" PRINT_INPUT\n") + f.write(" PRINT_FLOWS\n") + f.write(" OBS6 FILEIN uzf01.ghb.obs\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" bhead NETCDF\n") + f.write(" cond NETCDF\n") + f.write("END period 1\n") + +# TODO need extended modflow 6 to run this simulation +# run the netcdf sim +# success, buff = sim.run_simulation(silent=True, report=True) +# assert success, pformat(buff) diff --git a/.docs/Notebooks/netcdf02_tutorial.py b/.docs/Notebooks/netcdf02_tutorial.py new file mode 100644 index 000000000..4df8cbd29 --- /dev/null +++ b/.docs/Notebooks/netcdf02_tutorial.py @@ -0,0 +1,427 @@ +import sys +from pathlib import Path +from pprint import pformat +from tempfile import TemporaryDirectory + +import numpy as np +import xarray as xr + +import flopy + +print(sys.version) +print(f"flopy version: {flopy.__version__}") +# - + +DNODATA = 3.0e30 + + +def create_sim(ws): + name = "flow" + gwfname = name + sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name="mf6") + tdis_rc = [(100.0, 1, 1.0), (100.0, 1, 1.0)] + nper = len(tdis_rc) + tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc) + + gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True) + + # ims + hclose = 1.0e-6 + rclose = 1.0e-6 + 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 = 1 + nrow = 10 + ncol = 10 + delr = 10.0 + delc = 10.0 + top = 100.0 + botm = 0.0 + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + ic = flopy.mf6.ModflowGwfic(gwf, strt=100.0) + + npf = flopy.mf6.ModflowGwfnpf( + gwf, + xt3doptions=False, + save_flows=True, + save_specific_discharge=True, + save_saturation=True, + icelltype=[1], + k=10.0, + ) + + sto_on = False + if sto_on: + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=[1], + ss=1.0e-5, + sy=0.3, + steady_state={0: True}, + transient={0: 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 = False + if rch_on: + rch = flopy.mf6.ModflowGwfrcha(gwf, recharge={0: 4.79e-3}, pname="RCH-1") + + # wel + q = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + welconc = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + for i in range(nrow): + q[0, i, 0] = 100.0 + welconc[0, i, 0] = 100.0 + wel = flopy.mf6.ModflowGwfwelg( + gwf, + auxiliary=["concentration"], + pname="WEL-1", + q=q, + aux=welconc, + ) + + # ghb + rows = [0, 1, 2, 3] + for ipak, i in enumerate(rows): + fname = f"flow.{ipak + 1}.ghb" + pname = f"GHB-{ipak + 1}" + bhead = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + cond = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + conc = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + bhead[0, i, ncol - 1] = 50.0 + cond[0, i, ncol - 1] = 1000.0 + conc[0, i, ncol - 1] = 100.0 + flopy.mf6.ModflowGwfghbg( + gwf, + auxiliary=["concentration"], + filename=f"{fname}g", + pname=pname, + bhead=bhead, + cond=cond, + aux=conc, + ) + + # riv + rows = [4, 5, 6] + for ipak, i in enumerate(rows): + fname = f"flow.{ipak + 1}.riv" + pname = f"RIV-{ipak + 1}" + stage = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + cond = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + rbot = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + conc = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + stage[0, i, ncol - 1] = 50.0 + cond[0, i, ncol - 1] = 1000.0 + rbot[0, i, ncol - 1] = 0.0 + conc[0, i, ncol - 1] = 100.0 + riv = flopy.mf6.ModflowGwfrivg( + gwf, + auxiliary=["concentration"], + filename=f"{fname}g", + pname=pname, + stage=stage, + cond=cond, + rbot=rbot, + aux=[conc], + ) + + # drn + rows = [7, 8, 9] + for ipak, i in enumerate(rows): + fname = f"flow.{ipak + 1}.drn" + pname = f"DRN-{ipak + 1}" + elev = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + cond = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + conc = np.full((nlay, nrow, ncol), DNODATA, dtype=float) + elev[0, i, ncol - 1] = 50.0 + cond[0, i, ncol - 1] = 1000.0 + conc[0, i, ncol - 1] = 100.0 + drn = flopy.mf6.modflow.ModflowGwfdrng( + gwf, + auxiliary=["concentration"], + filename=f"{fname}g", + pname=pname, + elev=elev, + cond=cond, + aux=[conc], + ) + + return sim + + +def add_netcdf_vars(dataset, nc_info, dimmap): + def _data_shape(shape): + dims_l = [] + for d in shape: + dims_l.append(dimmap[d]) + + return dims_l + + for v in nc_info: + varname = nc_info[v]["varname"] + data = np.full( + _data_shape(nc_info[v]["nc_shape"]), + nc_info[v]["attrs"]["_FillValue"], + dtype=nc_info[v]["nc_type"], + ) + var_d = {varname: (nc_info[v]["nc_shape"], data)} + dataset = dataset.assign(var_d) + for a in nc_info[v]["attrs"]: + dataset[varname].attrs[a] = nc_info[v]["attrs"][a] + + return dataset + + +temp_dir = TemporaryDirectory() +workspace = Path(temp_dir.name) + +# run the non-netcdf simulation +sim = create_sim(workspace) +sim.write_simulation() +success, buff = sim.run_simulation(silent=True, report=True) +assert success, pformat(buff) + +# create directory for netcdf sim +sim.set_sim_path(workspace / "netcdf") +gwf = sim.get_model("flow") +gwf.name_file.nc_filerecord = "flow.structured.nc" +sim.write_simulation() + +# create the netcdf dataset +ds = xr.Dataset() + +# get model netcdf info +nc_info = gwf.netcdf_info() + +# update dataset with required attributes +for a in nc_info["attrs"]: + ds.attrs[a] = nc_info["attrs"][a] + +# get dim info from modelgrid +dis = gwf.modelgrid +xoff = dis.xoffset +yoff = dis.yoffset +x = xoff + dis.xycenters[0] +y = yoff + dis.xycenters[1] +z = [float(x) for x in range(1, dis.nlay + 1)] +nstp = sum(gwf.modeltime.nstp) +time = gwf.modeltime.tslen +nlay = dis.nlay +nrow = dis.nrow +ncol = dis.ncol +dimmap = {"time": nstp, "z": nlay, "y": nrow, "x": ncol} + +# create coordinate vars +var_d = {"time": (["time"], time), "z": (["z"], z), "y": (["y"], y), "x": (["x"], x)} +ds = ds.assign(var_d) + +# shape list for data arrays +shape = ["time", "z", "y", "x"] + +# update for welg +welg = gwf.get_package("wel-1") +nc_info = welg.netcdf_info() + +# create welg dataset variables +ds = add_netcdf_vars(ds, nc_info, dimmap) + +# update q netcdf array from flopy perioddata +for p in welg.q.get_data(): + if welg.q.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds["wel-1_q"].values[istp] = welg.q.get_data()[p] + +# update conc netcdf array from flopy perioddata +for p in welg.aux.get_data(): + if welg.aux.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds["wel-1_concentration"].values[istp] = welg.aux.get_data()[p][0] + +# update welg input to read from netcdf +with open(workspace / "netcdf" / "flow.welg", "w") as f: + f.write("BEGIN options\n") + f.write(" READARRAYGRID\n") + f.write(" auxiliary CONCENTRATION\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" q NETCDF\n") + f.write(" concentration NETCDF\n") + f.write("END period 1\n\n") + +# update for ghbg +for n in range(4): + ip = n + 1 + + # get ghbg package netcdf info + ghbg = gwf.get_package(f"ghb-{ip}") + nc_info = ghbg.netcdf_info() + + # create ghbg dataset variables + ds = add_netcdf_vars(ds, nc_info, dimmap) + + # update bhead netcdf array from flopy perioddata + for p in ghbg.bhead.get_data(): + if ghbg.bhead.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"ghb-{ip}_bhead"].values[istp] = ghbg.bhead.get_data()[p] + + # update cond netcdf array from flopy perioddata + for p in ghbg.cond.get_data(): + if ghbg.cond.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"ghb-{ip}_cond"].values[istp] = ghbg.cond.get_data()[p] + + # update conc netcdf array from flopy perioddata + for p in ghbg.aux.get_data(): + if ghbg.aux.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"ghb-{ip}_concentration"].values[istp] = ghbg.aux.get_data()[p][0] + + # update ghbg input to read from netcdf + with open(workspace / "netcdf" / f"flow.{ip}.ghbg", "w") as f: + f.write("BEGIN options\n") + f.write(" READARRAYGRID\n") + f.write(" auxiliary CONCENTRATION\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" bhead NETCDF\n") + f.write(" cond NETCDF\n") + f.write(" concentration NETCDF\n") + f.write("END period 1\n\n") + + +# update for rivg +for n in range(3): + ip = n + 1 + + # get rivg package netcdf info + rivg = gwf.get_package(f"riv-{ip}") + nc_info = rivg.netcdf_info() + + # create rivg dataset variables + ds = add_netcdf_vars(ds, nc_info, dimmap) + + # update stage netcdf array from flopy perioddata + for p in rivg.stage.get_data(): + if rivg.stage.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"riv-{ip}_stage"].values[istp] = rivg.stage.get_data()[p] + + # update cond netcdf array from flopy perioddata + for p in rivg.cond.get_data(): + if rivg.cond.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"riv-{ip}_cond"].values[istp] = rivg.cond.get_data()[p] + + # update rbot netcdf array from flopy perioddata + for p in rivg.rbot.get_data(): + if rivg.rbot.get_data()[p] is not None: + ds[f"riv-{ip}_rbot"].values[p] = rivg.rbot.get_data()[p] + + # update conc netcdf array from flopy perioddata + for p in rivg.aux.get_data(): + if rivg.aux.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"riv-{ip}_concentration"].values[istp] = rivg.aux.get_data()[p][0] + + # update rivg input to read from netcdf + with open(workspace / "netcdf" / f"flow.{ip}.rivg", "w") as f: + f.write("BEGIN options\n") + f.write(" READARRAYGRID\n") + f.write(" auxiliary CONCENTRATION\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" stage NETCDF\n") + f.write(" cond NETCDF\n") + f.write(" rbot NETCDF\n") + f.write(" concentration NETCDF\n") + f.write("END period 1\n\n") + + +# update for drng +for n in range(3): + ip = n + 1 + + # get drng package netcdf info + drng = gwf.get_package(f"drn-{ip}") + nc_info = drng.netcdf_info() + + # create drng dataset variables + ds = add_netcdf_vars(ds, nc_info, dimmap) + + # update elev netcdf array from flopy perioddata + for p in drng.elev.get_data(): + if drng.elev.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"drn-{ip}_elev"].values[istp] = drng.elev.get_data()[p] + + # update cond netcdf array from flopy perioddata + for p in drng.cond.get_data(): + if drng.cond.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"drn-{ip}_cond"].values[istp] = drng.cond.get_data()[p] + + # update conc netcdf array from flopy perioddata + for p in drng.aux.get_data(): + if drng.aux.get_data()[p] is not None: + istp = sum(gwf.modeltime.nstp[0:p]) + ds[f"drn-{ip}_concentration"].values[istp] = drng.aux.get_data()[p][0] + + # update drng input to read from netcdf + with open(workspace / "netcdf" / f"flow.{ip}.drng", "w") as f: + f.write("BEGIN options\n") + f.write(" READARRAYGRID\n") + f.write(" auxiliary CONCENTRATION\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" elev NETCDF\n") + f.write(" cond NETCDF\n") + f.write(" concentration NETCDF\n") + f.write("END period 1\n\n") + +# write the netcdf +ds.to_netcdf( + workspace / "netcdf/flow.structured.nc", format="NETCDF4", engine="netcdf4" +) + +# TODO need extended modflow 6 to run this simulation +# run the netcdf sim +# success, buff = sim.run_simulation(silent=True, report=True) +# assert success, pformat(buff) diff --git a/.github/workflows/rtd.yml b/.github/workflows/rtd.yml index e258d00b3..809201c51 100644 --- a/.github/workflows/rtd.yml +++ b/.github/workflows/rtd.yml @@ -79,6 +79,12 @@ jobs: echo $GITHUB_REF echo $GITHUB_EVENT_NAME + - name: Checkout MODFLOW 6 + uses: actions/checkout@v4 + with: + repository: MODFLOW-ORG/modflow6 + path: modflow6 + - name: Setup Micromamba uses: mamba-org/setup-micromamba@v1 with: @@ -135,6 +141,10 @@ jobs: ostag: mac subset: triangle + - name: Update FloPy packages + working-directory: modflow6 + run: python -m flopy.mf6.utils.generate_classes --dfnpath doc/mf6io/mf6ivar/dfn + - name: Run tutorial and example notebooks working-directory: autotest run: | diff --git a/etc/environment.yml b/etc/environment.yml index bcb8a08e3..d5dc32b5d 100644 --- a/etc/environment.yml +++ b/etc/environment.yml @@ -16,8 +16,7 @@ dependencies: - Jinja2>=3.0 - tomli - tomli-w - - pip: - - git+https://github.com/MODFLOW-ORG/modflow-devtools.git + - modflow-devtools==1.7.0 # lint - cffconvert @@ -31,8 +30,8 @@ dependencies: - jupyter - jupyter_client>=8.4.0 - jupytext + - modflow-devtools==1.7.0 - pip: - - git+https://github.com/MODFLOW-ORG/modflow-devtools.git - tach - pytest!=8.1.0 - pytest-benchmark @@ -64,3 +63,4 @@ dependencies: - vtk - xmipy - h5py + - xarray diff --git a/flopy/mf6/data/mfstructure.py b/flopy/mf6/data/mfstructure.py index 07115d12d..2d9ff246d 100644 --- a/flopy/mf6/data/mfstructure.py +++ b/flopy/mf6/data/mfstructure.py @@ -594,6 +594,7 @@ def __init__(self): self.parameter_name = None self.one_per_pkg = False self.jagged_array = None + self.netcdf = False def set_value(self, line, common): arr_line = line.strip().split() @@ -771,6 +772,8 @@ def set_value(self, line, common): self.one_per_pkg = bool(arr_line[1]) elif arr_line[0] == "jagged_array": self.jagged_array = arr_line[1] + elif arr_line[0] == "netcdf": + self.netcdf = arr_line[1] def get_type_string(self): return f"[{self.type_string}]" @@ -1089,6 +1092,7 @@ def __init__(self, data_item, model_data, package_type, dfn_list): or "nodes" in data_item.shape or len(data_item.layer_dims) > 1 ) + self.netcdf = data_item.netcdf self.num_data_items = len(data_item.data_items) self.record_within_record = False self.file_data = False diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index 0abc25cbd..c190522fc 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -2203,3 +2203,49 @@ def _resolve_idomain(idomain, botm): else: return np.ones_like(botm) return idomain + + @staticmethod + def netcdf_attrs(mname, mtype, grid_type, mesh=None): + """Return dictionary of dataset (model) scoped attributes + Parameters + ---------- + mname : str + model name + mtype : str + model type + grid_type: + DiscretizationType + mesh : str + mesh type if dataset is ugrid compliant + """ + attrs = { + "modflow_grid": "", + "modflow_model": "", + } + if grid_type == DiscretizationType.DIS: + attrs["modflow_grid"] = "STRUCTURED" + elif grid_type == DiscretizationType.DISV: + attrs["modflow_grid"] = "VERTEX" + + attrs["modflow_model"] = f"{mname.upper()}: MODFLOW 6 {mtype.upper()} model" + + # supported => LAYERED + if mesh: + attrs["mesh"] = mesh + + return attrs + + def netcdf_info(self, mesh=None): + """Return dictionary of dataset (model) scoped attributes + Parameters + ---------- + mesh : str + mesh type if dataset is ugrid compliant + """ + attrs = MFModel.netcdf_attrs( + self.name, self.model_type, self.get_grid_type(), mesh + ) + + res_d = {} + res_d["attrs"] = attrs + return res_d diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index 4ee9d6175..843931800 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -3422,6 +3422,162 @@ def plot(self, **kwargs): axes = PlotUtilities._plot_package_helper(self, **kwargs) return axes + @staticmethod + def _add_netcdf_entries( + attrs, mname, pname, data_item, auxiliary=None, mesh=None, nlay=1 + ): + DNODATA = 3.0e30 # MF6 DNODATA constant + FILLNA_INT32 = np.int32(-2147483647) # netcdf-fortran NF90_FILL_INT + FILLNA_DBL = 9.96920996838687e36 # netcdf-fortran NF90_FILL_DOUBLE + + if auxiliary: + auxnames = auxiliary + else: + auxnames = [] + + def _add_entry(tagname, iaux=None, layer=None): + + # netcdf variable dictionary + a = {} + + # set dict key and netcdf variable name + key = tagname + name = f"{pname}" + if iaux is not None: + key = f"{key}/{iaux}" + name = f"{name}_{auxiliary[iaux]}" + else: + name = f"{name}_{tagname}" + if layer is not None: + key = f"{key}/layer{layer}" + name = f"{name}_l{layer}" + + # add non-attrs to dictionary + a["varname"] = name.lower() + if (data_item.type) == DatumType.integer: + a["nc_type"] = np.int32 + elif (data_item.type) == DatumType.double_precision: + a["nc_type"] = np.float64 + dims = [] + if data_item.shape[0] == 'nodes': + if data_item.block_name == "griddata": + dims += ["x", "y", "z"] + elif data_item.block_name == "period": + dims += ["x", "y", "z", "time"] + else: + dimmap = {"nlay": "z", "nrow": "y", "ncol": "x"} + for s in data_item.shape: + for k, v in dimmap.items(): + s = s.replace(k, v) + dims.append(s) + a["nc_shape"] = dims[::-1] + + # add variable attributes dictionary + a["attrs"] = {} + a["attrs"]["modflow_input"] = (f"{mname}/{pname}/{tagname}").upper() + if iaux is not None: + a["attrs"]["modflow_iaux"] = iaux + 1 + if layer is not None: + a["attrs"]["layer"] = layer + if (data_item.type) == DatumType.integer: + a["attrs"]["_FillValue"] = FILLNA_INT32 + elif (data_item.type) == DatumType.double_precision: + if data_item.block_name == "griddata": + a["attrs"]["_FillValue"] = FILLNA_DBL + elif data_item.block_name == "period": + a["attrs"]["_FillValue"] = DNODATA + + # set dictionary + attrs[key] = a + + if data_item.layered and mesh == "LAYERED": + if data_item.name == "aux" or data_item.name == "auxvar": + for n, auxname in enumerate(auxnames): + for l in range(nlay): + _add_entry(data_item.name, n, l + 1) + else: + for l in range(nlay): + _add_entry(data_item.name, layer=l + 1) + else: + if data_item.name == "aux" or data_item.name == "auxvar": + for n, auxname in enumerate(auxnames): + _add_entry(data_item.name, iaux=n) + else: + _add_entry(data_item.name) + + @staticmethod + def netcdf_attrs(mtype, ptype, auxiliary=None, mesh=None, nlay=1): + from .data.mfstructure import DfnPackage, MFSimulationStructure + + attrs = {} + sim_struct = MFSimulationStructure() + + for package in MFPackage.__subclasses__(): + sim_struct.process_dfn(DfnPackage(package)) + p = DfnPackage(package) + c, sc = p.dfn_file_name.split(".")[0].split("-") + if c == mtype.lower() and sc == ptype.lower(): + sim_struct.add_package(p, model_file=False) + exit + + if ptype.lower() in sim_struct.package_struct_objs: + pso = sim_struct.package_struct_objs[ptype.lower()] + if pso.multi_package_support: + pname = f"<{ptype}name>" + else: + pname = ptype + for key, block in pso.blocks.items(): + if key != "griddata" and key != "period": + continue + for d in block.data_structures: + if block.data_structures[d].netcdf: + MFPackage._add_netcdf_entries( + attrs, + f"<{mtype}name>", + pname, + block.data_structures[d], + auxiliary, + mesh, + nlay, + ) + + res_d = {} + for k in list(attrs): + res_d[k] = attrs[k]["attrs"] + + return res_d + + def netcdf_info(self, mesh=None): + attrs = {} + + if self.dimensions.get_aux_variables(): + auxnames = list(self.dimensions.get_aux_variables()[0]) + if len(auxnames) and auxnames[0] == "auxiliary": + auxnames.pop(0) + else: + auxnames = [] + + for key, block in self.blocks.items(): + if key != "griddata" and key != "period": + continue + for dataset in block.datasets.values(): + if isinstance(dataset, mfdataarray.MFArray): + for index, data_item in enumerate( + dataset.structure.data_item_structures + ): + if dataset.structure.netcdf and dataset.has_data(): + MFPackage._add_netcdf_entries( + attrs, + self.model_name, + self.package_name, + dataset.structure, + auxnames, + mesh, + self.model_or_sim.modelgrid.nlay, + ) + + return attrs + class MFChildPackages: """