|
| 1 | +"""test_gwf_maw_pkgdata_ts.py |
| 2 | +
|
| 3 | +Focused tests for MAW PACKAGEDATA IDM enabled time-series linkage for two |
| 4 | +distinct TS mechanisms: |
| 5 | + Case 0: STRT package data timeseries variable |
| 6 | + Case 1: AUX package data timeseries variable |
| 7 | +
|
| 8 | +Model geometry: 1-layer, 1-row, 3-column; CHD at flanking columns; single |
| 9 | + MAW well in the centre column; 3 steady-state stress periods, 1 step each. |
| 10 | +""" |
| 11 | + |
| 12 | +import os |
| 13 | + |
| 14 | +import flopy |
| 15 | +import numpy as np |
| 16 | +import pytest |
| 17 | +from flopy.utils.compare import eval_bud_diff |
| 18 | +from framework import TestFramework |
| 19 | + |
| 20 | +paktest = "maw" |
| 21 | +cases = ["maw_strt_pd_ts", "maw_aux_pd_ts"] |
| 22 | + |
| 23 | +# Shared geometry and temporal discretisation |
| 24 | +nper = 3 |
| 25 | +perlen = 1.0 |
| 26 | +period_data = [(perlen, 1, 1.0)] * nper # steady-state, 1 timestep each |
| 27 | + |
| 28 | +nlay, nrow, ncol = 1, 1, 3 |
| 29 | +delr = delc = 100.0 |
| 30 | +top = 200.0 |
| 31 | +botm = [0.0] |
| 32 | +k11 = 1.0 |
| 33 | +chd_head = 100.0 |
| 34 | + |
| 35 | +# MAW geometry |
| 36 | +radius = 0.1 |
| 37 | +well_bot = -10.0 |
| 38 | + |
| 39 | +# Per-period values — chosen to be distinct so any failure to update is visible |
| 40 | +strt_vals = [90.0, 95.0, 85.0] |
| 41 | +conc_vals = [10.0, 20.0, 30.0] |
| 42 | + |
| 43 | +# TS times: entries at period-end times so LINEAREND gives exact per-period values. |
| 44 | +# tsmanager%ad() is called with totim = end of current timestep (1.0, 2.0, 3.0). |
| 45 | +ts_times = [0.0, 1.0, 2.0, 3.0] |
| 46 | +ts_strt = [strt_vals[0]] + strt_vals # [90, 90, 95, 85] |
| 47 | +ts_conc = [conc_vals[0]] + conc_vals # [10, 10, 20, 30] |
| 48 | + |
| 49 | + |
| 50 | +def _base_sim(ws, name): |
| 51 | + """Return a minimal GWF sim ready for a MAW package to be attached.""" |
| 52 | + sim = flopy.mf6.MFSimulation(sim_name=name, exe_name="mf6", sim_ws=ws) |
| 53 | + flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=period_data) |
| 54 | + flopy.mf6.ModflowIms(sim, print_option="NONE") |
| 55 | + gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) |
| 56 | + flopy.mf6.ModflowGwfdis( |
| 57 | + gwf, |
| 58 | + nlay=nlay, |
| 59 | + nrow=nrow, |
| 60 | + ncol=ncol, |
| 61 | + delr=delr, |
| 62 | + delc=delc, |
| 63 | + top=top, |
| 64 | + botm=botm, |
| 65 | + ) |
| 66 | + flopy.mf6.ModflowGwfnpf(gwf, k=k11) |
| 67 | + flopy.mf6.ModflowGwfic(gwf, strt=chd_head) |
| 68 | + flopy.mf6.ModflowGwfchd( |
| 69 | + gwf, |
| 70 | + stress_period_data=[[(0, 0, 0), chd_head], [(0, 0, 2), chd_head]], |
| 71 | + ) |
| 72 | + flopy.mf6.ModflowGwfoc( |
| 73 | + gwf, |
| 74 | + budget_filerecord=f"{name}.cbc", |
| 75 | + saverecord=[("BUDGET", "ALL")], |
| 76 | + ) |
| 77 | + return sim, gwf |
| 78 | + |
| 79 | + |
| 80 | +# Case 0 — STRT PACKAGEDATA TS |
| 81 | +def _get_strt_ts(ws, name): |
| 82 | + """STRT as PACKAGEDATA TS; STATUS CONSTANT; HEAD obs recorded.""" |
| 83 | + sim, gwf = _base_sim(ws, name) |
| 84 | + maw = flopy.mf6.ModflowGwfmaw( |
| 85 | + gwf, |
| 86 | + nmawwells=1, |
| 87 | + budget_filerecord=f"{name}.{paktest}.cbc", |
| 88 | + print_input=True, |
| 89 | + packagedata=[(0, radius, well_bot, "strt_ts", "THIEM", 1)], |
| 90 | + connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)], |
| 91 | + perioddata={0: [(0, "status", "constant")]}, |
| 92 | + pname="maw-1", |
| 93 | + ) |
| 94 | + maw.obs.initialize( |
| 95 | + filename=f"{name}.{paktest}.obs", |
| 96 | + digits=10, |
| 97 | + print_input=True, |
| 98 | + continuous={"maw_head.csv": [("well_head", "HEAD", (0,))]}, |
| 99 | + ) |
| 100 | + maw.ts.initialize( |
| 101 | + filename=f"{name}.{paktest}.ts", |
| 102 | + timeseries=list(zip(ts_times, ts_strt)), |
| 103 | + time_series_namerecord=["strt_ts"], |
| 104 | + interpolation_methodrecord=["linearend"], |
| 105 | + ) |
| 106 | + return sim |
| 107 | + |
| 108 | + |
| 109 | +# Case 1 — PACKAGEDATA AUX TS |
| 110 | +def _get_aux_period(ws, name): |
| 111 | + """Reference: AUX changes via PERIOD AUXILIARY each period; STRT fixed.""" |
| 112 | + sim, gwf = _base_sim(ws, name) |
| 113 | + maw_spd = { |
| 114 | + 0: [(0, "status", "constant"), (0, "AUXILIARY", "concentration", conc_vals[0])], |
| 115 | + 1: [(0, "AUXILIARY", "concentration", conc_vals[1])], |
| 116 | + 2: [(0, "AUXILIARY", "concentration", conc_vals[2])], |
| 117 | + } |
| 118 | + flopy.mf6.ModflowGwfmaw( |
| 119 | + gwf, |
| 120 | + nmawwells=1, |
| 121 | + auxiliary=["concentration"], |
| 122 | + budget_filerecord=f"{name}.{paktest}.cbc", |
| 123 | + packagedata=[(0, radius, well_bot, strt_vals[0], "THIEM", 1, conc_vals[0])], |
| 124 | + connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)], |
| 125 | + perioddata=maw_spd, |
| 126 | + pname="maw-1", |
| 127 | + ) |
| 128 | + return sim |
| 129 | + |
| 130 | + |
| 131 | +def _get_aux_pkgdata_ts(ws, name): |
| 132 | + """TS: AUX changes via PACKAGEDATA AUX TS; STRT fixed (same as reference).""" |
| 133 | + sim, gwf = _base_sim(ws, name) |
| 134 | + maw = flopy.mf6.ModflowGwfmaw( |
| 135 | + gwf, |
| 136 | + nmawwells=1, |
| 137 | + auxiliary=["concentration"], |
| 138 | + budget_filerecord=f"{name}.{paktest}.cbc", |
| 139 | + packagedata=[(0, radius, well_bot, strt_vals[0], "THIEM", 1, "conc_ts")], |
| 140 | + connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)], |
| 141 | + perioddata={0: [(0, "status", "constant")]}, |
| 142 | + pname="maw-1", |
| 143 | + ) |
| 144 | + maw.ts.initialize( |
| 145 | + filename=f"{name}.{paktest}.ts", |
| 146 | + timeseries=list(zip(ts_times, ts_conc)), |
| 147 | + time_series_namerecord=["conc_ts"], |
| 148 | + interpolation_methodrecord=["linearend"], |
| 149 | + ) |
| 150 | + return sim |
| 151 | + |
| 152 | + |
| 153 | +# Build the model |
| 154 | +def build_models(idx, test): |
| 155 | + name = cases[idx] |
| 156 | + ws0 = test.workspace |
| 157 | + ws1 = os.path.join(test.workspace, "mf6") |
| 158 | + if idx == 0: |
| 159 | + # Case 0: baseline is a fixed-STRT model (not compared); TS model |
| 160 | + # is in mf6/ sub-directory where HEAD obs are checked. |
| 161 | + return _get_strt_ts(ws0, name), _get_strt_ts(ws1, name) |
| 162 | + else: |
| 163 | + # Case 1: baseline uses PERIOD AUXILIARY; TS uses PACKAGEDATA AUX TS. |
| 164 | + return _get_aux_period(ws0, name), _get_aux_pkgdata_ts(ws1, name) |
| 165 | + |
| 166 | + |
| 167 | +# Check output |
| 168 | +def check_output(idx, test): |
| 169 | + name = os.path.basename(test.name) |
| 170 | + ws0 = test.workspace |
| 171 | + ws1 = os.path.join(test.workspace, "mf6") |
| 172 | + |
| 173 | + if idx == 0: |
| 174 | + # -- Case 0: assert that MAW HEAD obs (in mf6/) match STRT TS values -- |
| 175 | + obs_path = os.path.join(ws1, "maw_head.csv") |
| 176 | + obs = np.genfromtxt(obs_path, delimiter=",", names=True, skip_header=0) |
| 177 | + # One output time per stress period; values should equal strt_vals exactly |
| 178 | + # (LINEAREND at period-end totim gives the exact TS entry value) |
| 179 | + heads = obs["WELL_HEAD"] |
| 180 | + assert len(heads) == nper, f"expected {nper} obs records, got {len(heads)}" |
| 181 | + assert np.allclose(heads, strt_vals, atol=1e-6), ( |
| 182 | + f"STRT TS: well_head {heads} does not match expected {strt_vals}" |
| 183 | + ) |
| 184 | + else: |
| 185 | + # -- Case 1: GWF and MAW budgets must be identical between the two models -- |
| 186 | + ia = ( |
| 187 | + flopy.mf6.utils.MfGrdFile(os.path.join(ws0, f"{name}.dis.grb"))._datadict[ |
| 188 | + "IA" |
| 189 | + ] |
| 190 | + - 1 |
| 191 | + ) |
| 192 | + |
| 193 | + eval_bud_diff( |
| 194 | + os.path.join(ws0, f"{name}.cbc.cmp.out"), |
| 195 | + flopy.utils.CellBudgetFile( |
| 196 | + os.path.join(ws0, f"{name}.cbc"), precision="double" |
| 197 | + ), |
| 198 | + flopy.utils.CellBudgetFile( |
| 199 | + os.path.join(ws1, f"{name}.cbc"), precision="double" |
| 200 | + ), |
| 201 | + ia, |
| 202 | + ) |
| 203 | + eval_bud_diff( |
| 204 | + os.path.join(ws0, f"{name}.{paktest}.cbc.cmp.out"), |
| 205 | + flopy.utils.CellBudgetFile( |
| 206 | + os.path.join(ws0, f"{name}.{paktest}.cbc"), precision="double" |
| 207 | + ), |
| 208 | + flopy.utils.CellBudgetFile( |
| 209 | + os.path.join(ws1, f"{name}.{paktest}.cbc"), precision="double" |
| 210 | + ), |
| 211 | + ) |
| 212 | + |
| 213 | + |
| 214 | +@pytest.mark.parametrize("idx, name", enumerate(cases)) |
| 215 | +def test_mf6model(idx, name, function_tmpdir, targets): |
| 216 | + test = TestFramework( |
| 217 | + name=name, |
| 218 | + workspace=function_tmpdir, |
| 219 | + build=lambda t: build_models(idx, t), |
| 220 | + check=lambda t: check_output(idx, t), |
| 221 | + targets=targets, |
| 222 | + ) |
| 223 | + test.run() |
0 commit comments