Skip to content

Commit 0e14b2e

Browse files
committed
refactor advanced auxvar base integration with overrides
1 parent c764bf2 commit 0e14b2e

4 files changed

Lines changed: 272 additions & 48 deletions

File tree

autotest/test_gwf_maw_pkgdata_ts.py

Lines changed: 158 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@
44
distinct TS mechanisms:
55
Case 0: STRT package data timeseries variable
66
Case 1: AUX package data timeseries variable
7+
Case 2: verify PERIOD AUX data array override with PACKAGEDATA AUX TS
8+
Case 3: verify PERIOD AUX TS data array override with PACKAGEDATA AUX TS
79
810
Model geometry: 1-layer, 1-row, 3-column; CHD at flanking columns; single
911
MAW well in the centre column; 3 steady-state stress periods, 1 step each.
@@ -18,7 +20,7 @@
1820
from framework import TestFramework
1921

2022
paktest = "maw"
21-
cases = ["maw_strt_pd_ts", "maw_aux_pd_ts"]
23+
cases = ["maw_strt_pd_ts", "maw_aux_pd_ts", "maw_pd_oride", "maw_pd_ts_oride"]
2224

2325
# Shared geometry and temporal discretisation
2426
nper = 3
@@ -40,11 +42,16 @@
4042
strt_vals = [90.0, 95.0, 85.0]
4143
conc_vals = [10.0, 20.0, 30.0]
4244

45+
# PERIOD AUXILIARY override values — different from conc_vals so that a failure
46+
# to apply the PERIOD override is clearly distinguishable from the PKGDATA TS path.
47+
period_conc_vals = [40.0, 50.0, 60.0]
48+
4349
# TS times: entries at period-end times so LINEAREND gives exact per-period values.
4450
# tsmanager%ad() is called with totim = end of current timestep (1.0, 2.0, 3.0).
4551
ts_times = [0.0, 1.0, 2.0, 3.0]
4652
ts_strt = [strt_vals[0]] + strt_vals # [90, 90, 95, 85]
4753
ts_conc = [conc_vals[0]] + conc_vals # [10, 10, 20, 30]
54+
ts_period_conc = [period_conc_vals[0]] + period_conc_vals # [40, 40, 50, 60]
4855

4956

5057
def _base_sim(ws, name):
@@ -150,6 +157,143 @@ def _get_aux_pkgdata_ts(ws, name):
150157
return sim
151158

152159

160+
# Case 2 — PERIOD override PACKAGEDATA AUX TS
161+
# Reference model: PERIOD AUXILIARY with literal period_conc_vals; no PKGDATA TS.
162+
# TS model: PACKAGEDATA AUX TS (conc_vals) + same PERIOD AUXILIARY override.
163+
# If PERIOD correctly overrides, both models produce identical budget files.
164+
def _get_pd_ovrride_ref(ws, name):
165+
"""Reference: PERIOD AUXILIARY literal values; no PACKAGEDATA AUX TS."""
166+
sim, gwf = _base_sim(ws, name)
167+
maw_spd = {
168+
0: [
169+
(0, "status", "constant"),
170+
(0, "AUXILIARY", "concentration", period_conc_vals[0]),
171+
],
172+
1: [(0, "AUXILIARY", "concentration", period_conc_vals[1])],
173+
2: [(0, "AUXILIARY", "concentration", period_conc_vals[2])],
174+
}
175+
flopy.mf6.ModflowGwfmaw(
176+
gwf,
177+
nmawwells=1,
178+
auxiliary=["concentration"],
179+
budget_filerecord=f"{name}.{paktest}.cbc",
180+
packagedata=[
181+
(0, radius, well_bot, strt_vals[0], "THIEM", 1, period_conc_vals[0])
182+
],
183+
connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)],
184+
perioddata=maw_spd,
185+
pname="maw-1",
186+
)
187+
return sim
188+
189+
190+
def _get_pd_ovrride_ts(ws, name):
191+
"""TS: PKGDATA AUX TS (conc_vals) + PERIOD AUXILIARY override (period_conc_vals).
192+
193+
PERIOD must win: auxvar should equal period_conc_vals each period, not conc_vals.
194+
"""
195+
sim, gwf = _base_sim(ws, name)
196+
maw_spd = {
197+
0: [
198+
(0, "status", "constant"),
199+
(0, "AUXILIARY", "concentration", period_conc_vals[0]),
200+
],
201+
1: [(0, "AUXILIARY", "concentration", period_conc_vals[1])],
202+
2: [(0, "AUXILIARY", "concentration", period_conc_vals[2])],
203+
}
204+
maw = flopy.mf6.ModflowGwfmaw(
205+
gwf,
206+
nmawwells=1,
207+
auxiliary=["concentration"],
208+
budget_filerecord=f"{name}.{paktest}.cbc",
209+
packagedata=[(0, radius, well_bot, strt_vals[0], "THIEM", 1, "conc_ts")],
210+
connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)],
211+
perioddata=maw_spd,
212+
pname="maw-1",
213+
)
214+
maw.ts.initialize(
215+
filename=f"{name}.{paktest}.ts",
216+
timeseries=list(zip(ts_times, ts_conc)),
217+
time_series_namerecord=["conc_ts"],
218+
interpolation_methodrecord=["linearend"],
219+
)
220+
return sim
221+
222+
223+
# Case 3 — PERIOD AUX TS override PACKAGEDATA AUX TS
224+
# Reference model: PERIOD AUXILIARY TS (period_conc_ts) drives auxvar; no PKGDATA TS.
225+
# TS model: PACKAGEDATA AUX TS (conc_ts, conc_vals) +
226+
# PERIOD AUXILIARY TS (period_conc_ts, period_conc_vals).
227+
# If PERIOD TS correctly overrides, both models produce identical budget files.
228+
def _get_period_ts_ovrride_ref(ws, name):
229+
"""Reference: PERIOD AUXILIARY TS drives auxvar; no PACKAGEDATA AUX TS."""
230+
sim, gwf = _base_sim(ws, name)
231+
maw_spd = {
232+
0: [
233+
(0, "status", "constant"),
234+
(0, "AUXILIARY", "concentration", "period_conc_ts"),
235+
],
236+
1: [(0, "AUXILIARY", "concentration", "period_conc_ts")],
237+
2: [(0, "AUXILIARY", "concentration", "period_conc_ts")],
238+
}
239+
maw = flopy.mf6.ModflowGwfmaw(
240+
gwf,
241+
nmawwells=1,
242+
auxiliary=["concentration"],
243+
budget_filerecord=f"{name}.{paktest}.cbc",
244+
packagedata=[
245+
(0, radius, well_bot, strt_vals[0], "THIEM", 1, period_conc_vals[0])
246+
],
247+
connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)],
248+
perioddata=maw_spd,
249+
pname="maw-1",
250+
)
251+
maw.ts.initialize(
252+
filename=f"{name}.{paktest}.ts",
253+
timeseries=list(zip(ts_times, ts_period_conc)),
254+
time_series_namerecord=["period_conc_ts"],
255+
interpolation_methodrecord=["linearend"],
256+
)
257+
return sim
258+
259+
260+
def _get_period_ts_ovrride_ts(ws, name):
261+
"""TS: PKGDATA AUX TS (conc_ts) + PERIOD AUXILIARY TS (period_conc_ts).
262+
263+
PERIOD TS must win: auxvar should equal period_conc_vals, not conc_vals.
264+
Both TS series live in the same MAW TS file.
265+
"""
266+
sim, gwf = _base_sim(ws, name)
267+
maw_spd = {
268+
0: [
269+
(0, "status", "constant"),
270+
(0, "AUXILIARY", "concentration", "period_conc_ts"),
271+
],
272+
1: [(0, "AUXILIARY", "concentration", "period_conc_ts")],
273+
2: [(0, "AUXILIARY", "concentration", "period_conc_ts")],
274+
}
275+
maw = flopy.mf6.ModflowGwfmaw(
276+
gwf,
277+
nmawwells=1,
278+
auxiliary=["concentration"],
279+
budget_filerecord=f"{name}.{paktest}.cbc",
280+
packagedata=[(0, radius, well_bot, strt_vals[0], "THIEM", 1, "conc_ts")],
281+
connectiondata=[(0, 0, (0, 0, 1), top, well_bot, 0.0, 0.0)],
282+
perioddata=maw_spd,
283+
pname="maw-1",
284+
)
285+
# Both TS series in a single file: (time, conc_ts_val, period_conc_ts_val)
286+
maw.ts.initialize(
287+
filename=f"{name}.{paktest}.ts",
288+
timeseries=[
289+
(ts_times[i], ts_conc[i], ts_period_conc[i]) for i in range(len(ts_times))
290+
],
291+
time_series_namerecord=["conc_ts", "period_conc_ts"],
292+
interpolation_methodrecord=["linearend", "linearend"],
293+
)
294+
return sim
295+
296+
153297
# Build the model
154298
def build_models(idx, test):
155299
name = cases[idx]
@@ -159,9 +303,18 @@ def build_models(idx, test):
159303
# Case 0: baseline is a fixed-STRT model (not compared); TS model
160304
# is in mf6/ sub-directory where HEAD obs are checked.
161305
return _get_strt_ts(ws0, name), _get_strt_ts(ws1, name)
162-
else:
306+
elif idx == 1:
163307
# Case 1: baseline uses PERIOD AUXILIARY; TS uses PACKAGEDATA AUX TS.
164308
return _get_aux_period(ws0, name), _get_aux_pkgdata_ts(ws1, name)
309+
elif idx == 2:
310+
# Case 2: PERIOD override (literal) must win over PACKAGEDATA AUX TS.
311+
return _get_pd_ovrride_ref(ws0, name), _get_pd_ovrride_ts(ws1, name)
312+
else:
313+
# Case 3: PERIOD AUX TS must win over PACKAGEDATA AUX TS.
314+
return (
315+
_get_period_ts_ovrride_ref(ws0, name),
316+
_get_period_ts_ovrride_ts(ws1, name),
317+
)
165318

166319

167320
# Check output
@@ -182,7 +335,9 @@ def check_output(idx, test):
182335
f"STRT TS: well_head {heads} does not match expected {strt_vals}"
183336
)
184337
else:
185-
# -- Case 1: GWF and MAW budgets must be identical between the two models --
338+
# Cases 1-3: GWF and MAW budgets must be identical between reference and
339+
# TS model. For cases 2 and 3 this asserts that PERIOD (or PERIOD TS)
340+
# overrides PACKAGEDATA AUX TS.
186341
ia = (
187342
flopy.mf6.utils.MfGrdFile(os.path.join(ws0, f"{name}.dis.grb"))._datadict[
188343
"IA"

src/Model/GroundWaterFlow/gwf-maw.f90

Lines changed: 7 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,6 @@ module MawModule
5555
type :: MawInputContextType
5656
! -- packagedata
5757
real(DP), dimension(:), pointer, contiguous :: strt => null()
58-
integer(I4B), dimension(:), pointer, contiguous :: pkg_ifno => null()
5958
! -- period
6059
integer(I4B), pointer :: nbound => null()
6160
integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
@@ -773,7 +772,6 @@ end subroutine maw_input_init
773772
subroutine maw_input_destroy(this)
774773
class(MawInputContextType), intent(inout) :: this
775774
nullify (this%strt)
776-
nullify (this%pkg_ifno)
777775
nullify (this%nbound)
778776
nullify (this%ifno)
779777
nullify (this%setting)
@@ -824,7 +822,6 @@ subroutine maw_source_packagedata(this)
824822
call mem_setptr(radius, 'RADIUS', this%input_mempath)
825823
call mem_setptr(bottom, 'BOTTOM', this%input_mempath)
826824
call mem_setptr(this%input%strt, 'STRT', this%input_mempath)
827-
call mem_setptr(this%input%pkg_ifno, 'PACKAGEDATA_IFNO', this%input_mempath)
828825
call mem_setptr(condeqn, 'CONDEQN', this%input_mempath)
829826
call mem_setptr(ngwfnodes, 'NGWFNODES', this%input_mempath)
830827
if (this%inamedbound /= 0) then
@@ -1670,6 +1667,9 @@ subroutine maw_rp(this)
16701667
!
16711668
! -- check if input context has been updated for this stress period
16721669
if (this%iper == kper) then
1670+
!
1671+
! -- reset featureauxvar to pkg_auxvar baseline before period overrides
1672+
call this%set_auxvar_baseline()
16731673
!
16741674
! -- setup table to echo period data
16751675
if (this%iprpak /= 0) then
@@ -2019,18 +2019,14 @@ subroutine maw_ad(this)
20192019
integer(I4B) :: jj
20202020
integer(I4B) :: ibnd
20212021
integer(I4B) :: imaw
2022-
character(len=LINELENGTH) :: str
20232022
character(len=LENVARNAME) :: setting
20242023
!
2025-
! -- update AUX: bnd_ad() syncs featureauxvar from PACKAGEDATA AUX
2026-
! (timeseries supported) for every feature, establishing the
2027-
! per-timestep base value. The PERIOD loop below then overrides
2028-
! featureauxvar for any feature that has a PERIOD AUXILIARY entry,
2029-
! taking precedence over PACKAGEDATA.
2024+
! -- sync PACKAGEDATA AUX TS and advance observations;
2025+
! -- step auxvar period block overrides
20302026
call this%BndExtType%bnd_ad()
20312027
!
20322028
! -- sync PERIOD TS from input context
2033-
if (this%input%nbound > 0) then
2029+
if (this%ts_active .and. this%input%nbound > 0) then
20342030
do n = 1, this%input%nbound
20352031
imaw = this%input%ifno(n)
20362032
if (imaw < 1 .or. imaw > this%nmawwells) cycle
@@ -2044,17 +2040,6 @@ subroutine maw_ad(this)
20442040
this%well_head(imaw) = this%input%well_head(n)
20452041
this%xnewpak(imaw) = this%well_head(imaw)
20462042
end if
2047-
! AUXILIARY: overrides the PACKAGEDATA base set by bnd_ad() above
2048-
if (this%naux > 0) then
2049-
if (trim(setting) == 'PERIOD_AUXILIARY') then
2050-
str = this%input%auxname(n)
2051-
do jj = 1, this%naux
2052-
if (trim(str) /= trim(this%auxname(jj))) cycle
2053-
this%featureauxvar(jj, imaw) = this%input%auxval(n)
2054-
exit
2055-
end do
2056-
end if
2057-
end if
20582043
end do
20592044
end if
20602045
!
@@ -2076,7 +2061,7 @@ subroutine maw_ad(this)
20762061
! This unconditionally overwrites any WELL_HEAD set in the PERIOD loop
20772062
! above for CONSTANT wells.
20782063
do i = 1, this%nmawwells
2079-
n = this%input%pkg_ifno(i)
2064+
n = this%pkg_ifno(i)
20802065
this%xoldpak(n) = this%xnewpak(n)
20812066
this%xoldsto(n) = this%xsto(n)
20822067
if (this%iboundpak(n) < 0) then

0 commit comments

Comments
 (0)