Skip to content

Commit

Permalink
test: compare mp7 & prt results in notebook tests, misc cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
wpbonelli committed Aug 31, 2023
1 parent 26acaee commit e78ba7b
Show file tree
Hide file tree
Showing 9 changed files with 405 additions and 290 deletions.
51 changes: 0 additions & 51 deletions autotest/prt_test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,56 +54,6 @@ def check_track_data(
assert all(data_bin["ireason"] >= 0)


def to_mp7_format(data: Union[pd.DataFrame, np.recarray]) -> pd.DataFrame:
if isinstance(data, pd.DataFrame):
data = data.to_records(index=False)

mp7_dtypes = np.dtype(
[
("particleid", np.int32),
("particlegroup", np.int32),
("sequencenumber", np.int32),
("particleidloc", np.int32),
("time", np.float32),
("x", np.float32),
("y", np.float32),
("z", np.float32),
("k", np.int32),
("node", np.int32),
("xloc", np.float32),
("yloc", np.float32),
("zloc", np.float32),
("stressperiod", np.int32),
("timestep", np.int32),
]
)

return pd.DataFrame(
np.core.records.fromarrays(
[
data["irpt"],
data["iprp"],
np.zeros(
data.shape[0]
), # todo use sequencenumber passed explicitly as particle name
np.zeros(data.shape[0]),
data["t"],
data["x"],
data["y"],
data["z"],
data["ilay"], # todo add to PRT output?
data["icell"],
np.zeros(data.shape[0]),
np.zeros(data.shape[0]),
np.zeros(data.shape[0]),
data["kper"],
data["kstp"],
],
dtype=mp7_dtypes,
)
)


def check_budget_data(lst: os.PathLike, perlen=1, nper=1, nstp=1):
# load PRT model's list file
mflist = flopy.utils.mflistfile.ListBudget(
Expand Down Expand Up @@ -190,4 +140,3 @@ def get_partdata(grid, rpts):
timeoffset=0,
drape=0,
)

59 changes: 26 additions & 33 deletions autotest/test_prt_exg01.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,14 @@
import pytest
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.plot.plotutil import to_mp7_pathlines

from framework import TestFramework
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
has_default_boundnames,
to_mp7_format,
)
from simulation import TestSimulation

Expand Down Expand Up @@ -266,32 +266,31 @@ def test_mf6model(idx, name, function_tmpdir, targets):

# load mp7 pathline results
plf = PathlineFile(ws / mp7_pathline_file)
mp7_pldata = pd.DataFrame(
mp7_pls = pd.DataFrame(
plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
)
# convert zero-based to one-based
mp7_pldata["particleid"] = mp7_pldata["particleid"] + 1
mp7_pldata["particlegroup"] = mp7_pldata["particlegroup"] + 1
mp7_pldata["node"] = mp7_pldata["node"] + 1
mp7_pldata["k"] = mp7_pldata["k"] + 1
mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1
mp7_pls["node"] = mp7_pls["node"] + 1
mp7_pls["k"] = mp7_pls["k"] + 1

# load mf6 pathline results
mf6_pldata = pd.read_csv(ws / prt_track_csv_file).replace(
mf6_pls = pd.read_csv(ws / prt_track_csv_file).replace(
r"^\s*$", np.nan, regex=True
)

# make sure pathline dataframe has "name" column
assert "name" in mf6_pldata
assert "name" in mf6_pls

# check boundname values
if "bnms" in name:
# boundnames should be release point numbers (so pandas parses them as ints)
assert np.array_equal(
mf6_pldata["name"].to_numpy(), mf6_pldata["irpt"].to_numpy()
mf6_pls["name"].to_numpy(), mf6_pls["irpt"].to_numpy()
)
else:
# no boundnames given so check for defaults
assert pd.isna(mf6_pldata["name"]).all()
assert pd.isna(mf6_pls["name"]).all()

# check budget data were written to mf6 prt list file
check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
Expand Down Expand Up @@ -320,7 +319,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mf6_plines = mf6_pldata.groupby(["iprp", "irpt", "trelease"])
mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
pl.plot(
title="MF6 pathlines",
Expand All @@ -337,7 +336,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mp7_plines = mp7_pldata.groupby(["particleid"])
mp7_plines = mp7_pls.groupby(["particleid"])
for ipl, (pid, pl) in enumerate(mp7_plines):
pl.plot(
title="MP7 pathlines",
Expand All @@ -354,30 +353,24 @@ def test_mf6model(idx, name, function_tmpdir, targets):
plt.savefig(ws / f"test_{name}.png")

# convert mf6 pathlines to mp7 format
mf6_pldata_mp7 = to_mp7_format(mf6_pldata)
mf6_pls = to_mp7_pathlines(mf6_pls)

# sort both dataframes by particleid and time
mf6_pldata_mp7.sort_values(by=["particleid", "time"], inplace=True)
mp7_pldata.sort_values(by=["particleid", "time"], inplace=True)

# drop duplicate locations
# (mp7 includes a duplicate location at the end of each pathline??)
cols = ["particleid", "x", "y", "z", "time"]
mp7_pldata = mp7_pldata.drop_duplicates(subset=cols)
mf6_pldata_mp7 = mf6_pldata_mp7.drop_duplicates(subset=cols)
mf6_pls.sort_values(by=["particleid", "time"], inplace=True)
mp7_pls.sort_values(by=["particleid", "time"], inplace=True)

# drop columns for which there is no direct correspondence between mf6 and mp7
del mf6_pldata_mp7["sequencenumber"]
del mf6_pldata_mp7["particleidloc"]
del mf6_pldata_mp7["xloc"]
del mf6_pldata_mp7["yloc"]
del mf6_pldata_mp7["zloc"]
del mp7_pldata["sequencenumber"]
del mp7_pldata["particleidloc"]
del mp7_pldata["xloc"]
del mp7_pldata["yloc"]
del mp7_pldata["zloc"]
del mf6_pls["sequencenumber"]
del mf6_pls["particleidloc"]
del mf6_pls["xloc"]
del mf6_pls["yloc"]
del mf6_pls["zloc"]
del mp7_pls["sequencenumber"]
del mp7_pls["particleidloc"]
del mp7_pls["xloc"]
del mp7_pls["yloc"]
del mp7_pls["zloc"]

# compare mf6 / mp7 pathline data
assert mf6_pldata_mp7.shape == mp7_pldata.shape
assert np.allclose(mf6_pldata_mp7, mp7_pldata, atol=1e-3)
assert mf6_pls.shape == mp7_pls.shape
assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
55 changes: 27 additions & 28 deletions autotest/test_prt_fmi01.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,13 @@
import pytest
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.plot.plotutil import to_mp7_pathlines
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_partdata,
has_default_boundnames,
to_mp7_format,
)


Expand Down Expand Up @@ -339,25 +339,24 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):

# load mp7 pathline results
plf = PathlineFile(ws / mp7_pathline_file)
mp7_pldata = pd.DataFrame(
mp7_pls = pd.DataFrame(
plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
)
# convert zero-based to one-based indexing in mp7 results
mp7_pldata["particleid"] = mp7_pldata["particleid"] + 1
mp7_pldata["particlegroup"] = mp7_pldata["particlegroup"] + 1
mp7_pldata["node"] = mp7_pldata["node"] + 1
mp7_pldata["k"] = mp7_pldata["k"] + 1
mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1
mp7_pls["node"] = mp7_pls["node"] + 1
mp7_pls["k"] = mp7_pls["k"] + 1

# load mf6 pathline results
mf6_pldata = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
mf6_pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)

# make sure pathline df has "name" (boundname) column and default values
assert "name" in mf6_pldata
assert has_default_boundnames(mf6_pldata)
assert "name" in mf6_pls
assert has_default_boundnames(mf6_pls)

# make sure all mf6 pathline data have correct model and PRP index (1)
assert all_equal(mf6_pldata["imdl"], 1)
assert all_equal(mf6_pldata["iprp"], 1)
assert all_equal(mf6_pls["imdl"], 1)
assert all_equal(mf6_pls["iprp"], 1)

# check budget data were written to mf6 prt list file
check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
Expand Down Expand Up @@ -387,7 +386,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mf6_plines = mf6_pldata.groupby(["iprp", "irpt", "trelease"])
mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
pl.plot(
title="MF6 pathlines",
Expand All @@ -404,7 +403,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mp7_plines = mp7_pldata.groupby(["particleid"])
mp7_plines = mp7_pls.groupby(["particleid"])
for ipl, (pid, pl) in enumerate(mp7_plines):
pl.plot(
title="MP7 pathlines",
Expand All @@ -421,24 +420,24 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
plt.savefig(ws / f"test_{simname}.png")

# convert mf6 pathlines to mp7 format
mf6_pldata_mp7 = to_mp7_format(mf6_pldata)
mf6_pls = to_mp7_pathlines(mf6_pls)

# sort both dataframes by particleid and time
mf6_pldata_mp7.sort_values(by=["particleid", "time"], inplace=True)
mp7_pldata.sort_values(by=["particleid", "time"], inplace=True)
mf6_pls.sort_values(by=["particleid", "time"], inplace=True)
mp7_pls.sort_values(by=["particleid", "time"], inplace=True)

# drop columns for which there is no direct correspondence between mf6 and mp7
del mf6_pldata_mp7["sequencenumber"]
del mf6_pldata_mp7["particleidloc"]
del mf6_pldata_mp7["xloc"]
del mf6_pldata_mp7["yloc"]
del mf6_pldata_mp7["zloc"]
del mp7_pldata["sequencenumber"]
del mp7_pldata["particleidloc"]
del mp7_pldata["xloc"]
del mp7_pldata["yloc"]
del mp7_pldata["zloc"]
del mf6_pls["sequencenumber"]
del mf6_pls["particleidloc"]
del mf6_pls["xloc"]
del mf6_pls["yloc"]
del mf6_pls["zloc"]
del mp7_pls["sequencenumber"]
del mp7_pls["particleidloc"]
del mp7_pls["xloc"]
del mp7_pls["yloc"]
del mp7_pls["zloc"]

# compare mf6 / mp7 pathline data
assert mf6_pldata_mp7.shape == mp7_pldata.shape
assert np.allclose(mf6_pldata_mp7, mp7_pldata, atol=1e-3)
assert mf6_pls.shape == mp7_pls.shape
assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
Loading

0 comments on commit e78ba7b

Please sign in to comment.