|
| 1 | +""" |
| 2 | +Tests ability to transform release point coordinates from |
| 3 | +global, local, or local offset coordinate representations |
| 4 | +to model coordinates, as well as to transform output data |
| 5 | +from model to global cooordinates. |
| 6 | +
|
| 7 | +The grid is a 10x10 square with a single layer, |
| 8 | +the same flow system shown on the FloPy readme. |
| 9 | +
|
| 10 | +Test cases are defined for each coordinate system option. |
| 11 | +These are: |
| 12 | +
|
| 13 | + - local_xy |
| 14 | + - local_xy_offset |
| 15 | + - dev_global_xy |
| 16 | +
|
| 17 | +""" |
| 18 | + |
| 19 | +from pathlib import Path |
| 20 | + |
| 21 | +import flopy |
| 22 | +import matplotlib.cm as cm |
| 23 | +import matplotlib.pyplot as plt |
| 24 | +import numpy as np |
| 25 | +import pandas as pd |
| 26 | +import pytest |
| 27 | +from flopy.utils import PathlineFile |
| 28 | +from flopy.utils.binaryfile import HeadFile |
| 29 | +from flopy.utils.gridintersect import GridIntersect |
| 30 | +from shapely.geometry import Point |
| 31 | +from framework import TestFramework |
| 32 | +from prt_test_utils import ( |
| 33 | + FlopyReadmeCase, |
| 34 | + all_equal, |
| 35 | + check_budget_data, |
| 36 | + check_track_data, |
| 37 | + get_model_name, |
| 38 | + get_partdata, |
| 39 | + has_default_boundnames, |
| 40 | + DEFAULT_EXIT_SOLVE_TOL, |
| 41 | +) |
| 42 | + |
| 43 | +simname = "prtcrd" |
| 44 | +cases = [f"{simname}mdl", f"{simname}lcl", f"{simname}ofs"] # todo global |
| 45 | +nodes = list(range(100)) |
| 46 | + |
| 47 | + |
| 48 | +def get_start_points(): |
| 49 | + return np.add( |
| 50 | + np.transpose( |
| 51 | + np.array(np.meshgrid(range(10), range(10))).reshape(2, -1) |
| 52 | + ), |
| 53 | + 0.5, |
| 54 | + ) |
| 55 | + |
| 56 | + |
| 57 | +def get_prp(prt): |
| 58 | + gi = GridIntersect(prt.modelgrid, method="vertex", rtree=True) |
| 59 | + if "mdl" in prt.name: |
| 60 | + startpts = get_start_points() |
| 61 | + releasepts = [ |
| 62 | + (i, (0, *gi.intersect(Point(rpt))[0].cellids), *rpt, 0.5) |
| 63 | + for i, rpt in enumerate(startpts) |
| 64 | + ] |
| 65 | + elif "lcl" in prt.name: |
| 66 | + releasepts = [ |
| 67 | + (nn, tuple([*prt.modelgrid.get_lrc([nn])[0]]), 0.5, 0.5, 0.5) |
| 68 | + for nn in nodes |
| 69 | + ] |
| 70 | + elif "ofs" in prt.name: |
| 71 | + releasepts = [ |
| 72 | + (nn, tuple([*prt.modelgrid.get_lrc([nn])[0]]), 0.0, 0.0, 0.0) |
| 73 | + for nn in nodes |
| 74 | + ] |
| 75 | + elif "gbl" in prt.name: |
| 76 | + # todo global |
| 77 | + pass |
| 78 | + |
| 79 | + prp_track_file = f"{prt.name}_prt.prp.trk" |
| 80 | + prp_track_csv_file = f"{prt.name}_prt.prp.trk.csv" |
| 81 | + return flopy.mf6.ModflowPrtprp( |
| 82 | + prt, |
| 83 | + pname="prp1", |
| 84 | + filename=f"{prt.name}_1.prp", |
| 85 | + nreleasepts=len(releasepts), |
| 86 | + packagedata=releasepts, |
| 87 | + perioddata={0: ["FIRST"]}, |
| 88 | + track_filerecord=[prp_track_file], |
| 89 | + trackcsv_filerecord=[prp_track_csv_file], |
| 90 | + local_xy="lcl" in prt.name, |
| 91 | + local_xy_offset="ofs" in prt.name, |
| 92 | + exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, |
| 93 | + extend_tracking=True, |
| 94 | + ) |
| 95 | + |
| 96 | + |
| 97 | +def build_prt_sim(name, gwf_ws, prt_ws, mf6): |
| 98 | + # create simulation |
| 99 | + sim = flopy.mf6.MFSimulation( |
| 100 | + sim_name=name, |
| 101 | + exe_name=mf6, |
| 102 | + version="mf6", |
| 103 | + sim_ws=prt_ws, |
| 104 | + ) |
| 105 | + |
| 106 | + # create tdis package |
| 107 | + flopy.mf6.modflow.mftdis.ModflowTdis( |
| 108 | + sim, |
| 109 | + pname="tdis", |
| 110 | + time_units="DAYS", |
| 111 | + nper=FlopyReadmeCase.nper, |
| 112 | + perioddata=[ |
| 113 | + ( |
| 114 | + FlopyReadmeCase.perlen, |
| 115 | + FlopyReadmeCase.nstp, |
| 116 | + FlopyReadmeCase.tsmult, |
| 117 | + ) |
| 118 | + ], |
| 119 | + ) |
| 120 | + |
| 121 | + # create prt model |
| 122 | + prt_name = get_model_name(name, "prt") |
| 123 | + prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name, save_flows=True) |
| 124 | + |
| 125 | + # create prt discretization |
| 126 | + flopy.mf6.modflow.mfgwfdis.ModflowGwfdis( |
| 127 | + prt, |
| 128 | + pname="dis", |
| 129 | + nlay=FlopyReadmeCase.nlay, |
| 130 | + nrow=FlopyReadmeCase.nrow, |
| 131 | + ncol=FlopyReadmeCase.ncol, |
| 132 | + ) |
| 133 | + |
| 134 | + # create mip package |
| 135 | + flopy.mf6.ModflowPrtmip( |
| 136 | + prt, pname="mip", porosity=FlopyReadmeCase.porosity |
| 137 | + ) |
| 138 | + |
| 139 | + # create prp package |
| 140 | + prp = get_prp(prt) |
| 141 | + |
| 142 | + # create output control package |
| 143 | + prt_budget_file = f"{prt_name}.bud" |
| 144 | + prt_track_file = f"{prt_name}.trk" |
| 145 | + prt_track_csv_file = f"{prt_name}.trk.csv" |
| 146 | + flopy.mf6.ModflowPrtoc( |
| 147 | + prt, |
| 148 | + pname="oc", |
| 149 | + budget_filerecord=[prt_budget_file], |
| 150 | + track_filerecord=[prt_track_file], |
| 151 | + trackcsv_filerecord=[prt_track_csv_file], |
| 152 | + saverecord=[("BUDGET", "ALL")], |
| 153 | + ) |
| 154 | + |
| 155 | + # create the flow model interface |
| 156 | + gwf_name = get_model_name(name, "gwf") |
| 157 | + gwf_budget_file = gwf_ws / f"{gwf_name}.bud" |
| 158 | + gwf_head_file = gwf_ws / f"{gwf_name}.hds" |
| 159 | + flopy.mf6.ModflowPrtfmi( |
| 160 | + prt, |
| 161 | + packagedata=[ |
| 162 | + ("GWFHEAD", gwf_head_file), |
| 163 | + ("GWFBUDGET", gwf_budget_file), |
| 164 | + ], |
| 165 | + ) |
| 166 | + |
| 167 | + # add explicit model solution |
| 168 | + ems = flopy.mf6.ModflowEms( |
| 169 | + sim, |
| 170 | + pname="ems", |
| 171 | + filename=f"{prt_name}.ems", |
| 172 | + ) |
| 173 | + sim.register_solution_package(ems, [prt.name]) |
| 174 | + |
| 175 | + return sim |
| 176 | + |
| 177 | + |
| 178 | +def build_models(idx, test): |
| 179 | + gwf_ws = test.workspace / "gwf" |
| 180 | + prt_ws = test.workspace / "prt" |
| 181 | + gwf_sim = FlopyReadmeCase.get_gwf_sim( |
| 182 | + test.name, gwf_ws, test.targets["mf6"] |
| 183 | + ) |
| 184 | + prt_sim = build_prt_sim(test.name, gwf_ws, prt_ws, test.targets["mf6"]) |
| 185 | + return gwf_sim, prt_sim |
| 186 | + |
| 187 | + |
| 188 | +def check_output(idx, test): |
| 189 | + name = test.name |
| 190 | + gwf_ws = test.workspace |
| 191 | + prt_ws = test.workspace / "prt" |
| 192 | + gwf_sim, prt_sim = test.sims[:2] |
| 193 | + gwf = gwf_sim.get_model() |
| 194 | + prt = prt_sim.get_model() |
| 195 | + hds = gwf.output.head().get_data() |
| 196 | + mg = gwf.modelgrid |
| 197 | + pathlines = pd.read_csv(prt_ws / f"{prt.name}.trk.csv") |
| 198 | + startpts = pathlines[pathlines.ireason == 0] |
| 199 | + endpts = pathlines[pathlines.ireason == 3] |
| 200 | + |
| 201 | + # check start points |
| 202 | + assert np.allclose( |
| 203 | + np.sort(startpts[["x", "y"]].to_numpy(), axis=0), |
| 204 | + np.sort(np.array(get_start_points()), axis=0), |
| 205 | + ) |
| 206 | + |
| 207 | + # setup plot |
| 208 | + plot_data = False |
| 209 | + if plot_data: |
| 210 | + fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10)) |
| 211 | + for a in ax: |
| 212 | + a.set_aspect("equal") |
| 213 | + |
| 214 | + # plot pathline in map view |
| 215 | + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) |
| 216 | + pmv.plot_grid() |
| 217 | + pmv.plot_array(hds[0], alpha=0.1, cmap="jet") |
| 218 | + pmv.plot_pathline(pathlines) |
| 219 | + |
| 220 | + # view/save plot |
| 221 | + plt.show() |
| 222 | + plt.savefig(gwf_ws / f"test_{simname}.png") |
| 223 | + |
| 224 | + |
| 225 | +@pytest.mark.parametrize("idx, name", enumerate(cases)) |
| 226 | +def test_mf6model(idx, name, function_tmpdir, targets): |
| 227 | + test = TestFramework( |
| 228 | + name=name, |
| 229 | + workspace=function_tmpdir, |
| 230 | + build=lambda t: build_models(idx, t), |
| 231 | + check=lambda t: check_output(idx, t), |
| 232 | + targets=targets, |
| 233 | + compare=None, |
| 234 | + ) |
| 235 | + test.run() |
0 commit comments