Skip to content

Commit e5d19cc

Browse files
wpbonelliwpbonelli
committed
fix(prt): terminate particle if stop time precedes release time
1 parent 558c722 commit e5d19cc

2 files changed

Lines changed: 195 additions & 3 deletions

File tree

‎autotest/test_prt_stoptime_before_release.py‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
"""
2+
Regression test for https://github.com/MODFLOW-ORG/modflow6/issues/2941
3+
4+
If a PRP package's STOPTIME precedes a particle's release time, the particle
5+
is handed a negative tracking interval on release (tmax < ttrack), which used
6+
to segfault/crash (an FPE overflow trap in the analytic tracking formulas).
7+
MODFLOW 6 should instead warn and terminate the particle immediately without
8+
tracking it.
9+
10+
Two release configurations are tested, matching the two models attached to
11+
the issue:
12+
- relt: an explicit RELEASETIMES entry falls after STOPTIME
13+
- prd: a PERIOD block release (FIRST) fills forward into a later stress
14+
period whose start time falls after STOPTIME
15+
"""
16+
17+
import flopy
18+
import pytest
19+
from framework import TestFramework
20+
from prt_test_utils import FlopyReadmeCase, get_model_name
21+
22+
simname = "prtstpb4"
23+
cases = [f"{simname}rt", f"{simname}pd"]
24+
25+
# release point matching cell (0, 0, 0) in the FlopyReadmeCase grid
26+
releasepts = [[0, 0, 0, 0, 0.1, 9.1, 0.5]]
27+
28+
29+
def build_prt_sim(name, gwf_ws, prt_ws, mf6):
30+
# create simulation
31+
sim = flopy.mf6.MFSimulation(
32+
sim_name=name,
33+
exe_name=mf6,
34+
version="mf6",
35+
sim_ws=prt_ws,
36+
)
37+
38+
# two stress periods so a later-period release time/step can fall
39+
# after STOPTIME, which is set to expire during the first period
40+
flopy.mf6.modflow.mftdis.ModflowTdis(
41+
sim,
42+
pname="tdis",
43+
time_units="DAYS",
44+
nper=2,
45+
perioddata=[
46+
(FlopyReadmeCase.perlen, FlopyReadmeCase.nstp, FlopyReadmeCase.tsmult),
47+
(FlopyReadmeCase.perlen, FlopyReadmeCase.nstp, FlopyReadmeCase.tsmult),
48+
],
49+
)
50+
51+
# create prt model
52+
prt_name = get_model_name(name, "prt")
53+
prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
54+
55+
# create prt discretization
56+
flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
57+
prt,
58+
pname="dis",
59+
nlay=FlopyReadmeCase.nlay,
60+
nrow=FlopyReadmeCase.nrow,
61+
ncol=FlopyReadmeCase.ncol,
62+
top=FlopyReadmeCase.top,
63+
botm=FlopyReadmeCase.botm,
64+
)
65+
66+
# create mip package
67+
flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=FlopyReadmeCase.porosity)
68+
69+
# create prp package
70+
prp_track_file = f"{prt_name}.prp.trk"
71+
prp_track_csv_file = f"{prt_name}.prp.trk.csv"
72+
73+
if name.endswith("rt"):
74+
# release time (1.5) falls in period 2, after stoptime (0.5)
75+
perioddata = None
76+
nreleasetimes = 1
77+
releasetimes = [(1.5,)]
78+
else:
79+
# FIRST in period 1 releases at t=0.0 (before stoptime), but fills
80+
# forward into period 2, releasing again at t=1.0 (after stoptime)
81+
perioddata = {0: [("FIRST",)]}
82+
nreleasetimes = None
83+
releasetimes = None
84+
85+
flopy.mf6.ModflowPrtprp(
86+
prt,
87+
pname="prp1",
88+
filename=f"{prt_name}_1.prp",
89+
nreleasepts=len(releasepts),
90+
packagedata=releasepts,
91+
perioddata=perioddata,
92+
nreleasetimes=nreleasetimes,
93+
releasetimes=releasetimes,
94+
stoptime=0.5,
95+
track_filerecord=[prp_track_file],
96+
trackcsv_filerecord=[prp_track_csv_file],
97+
print_input=True,
98+
extend_tracking=True,
99+
)
100+
101+
# create output control package
102+
prt_track_file = f"{prt_name}.trk"
103+
prt_track_csv_file = f"{prt_name}.trk.csv"
104+
flopy.mf6.ModflowPrtoc(
105+
prt,
106+
pname="oc",
107+
track_filerecord=[prt_track_file],
108+
trackcsv_filerecord=[prt_track_csv_file],
109+
)
110+
111+
# create the flow model interface
112+
gwf_name = get_model_name(name, "gwf")
113+
gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
114+
gwf_head_file = gwf_ws / f"{gwf_name}.hds"
115+
flopy.mf6.ModflowPrtfmi(
116+
prt,
117+
packagedata=[
118+
("GWFHEAD", gwf_head_file),
119+
("GWFBUDGET", gwf_budget_file),
120+
],
121+
)
122+
123+
# add explicit model solution
124+
ems = flopy.mf6.ModflowEms(
125+
sim,
126+
pname="ems",
127+
filename=f"{prt_name}.ems",
128+
)
129+
sim.register_solution_package(ems, [prt.name])
130+
131+
return sim
132+
133+
134+
def build_models(test):
135+
gwf_sim = FlopyReadmeCase.get_gwf_sim(
136+
test.name, test.workspace, test.targets["mf6"]
137+
)
138+
# GWF sim also needs 2 stress periods to match the PRT model's TDIS
139+
tdis = gwf_sim.get_package("tdis")
140+
tdis.nper = 2
141+
tdis.perioddata = [
142+
(FlopyReadmeCase.perlen, FlopyReadmeCase.nstp, FlopyReadmeCase.tsmult),
143+
(FlopyReadmeCase.perlen, FlopyReadmeCase.nstp, FlopyReadmeCase.tsmult),
144+
]
145+
prt_sim = build_prt_sim(
146+
test.name,
147+
test.workspace,
148+
test.workspace / "prt",
149+
test.targets["mf6"],
150+
)
151+
return gwf_sim, prt_sim
152+
153+
154+
def check_output(test):
155+
# simulation should complete normally (not crash) despite stoptime
156+
# preceding the particle's release time, and should warn about it.
157+
# normalize whitespace since the message may be line-wrapped in the
158+
# list file.
159+
lst = " ".join((test.workspace / "prt" / "mfsim.lst").read_text().split())
160+
assert "already passed; terminating without tracking" in lst
161+
162+
163+
@pytest.mark.parametrize("name", cases)
164+
def test_mf6model(name, function_tmpdir, targets):
165+
test = TestFramework(
166+
name=name,
167+
workspace=function_tmpdir,
168+
build=build_models,
169+
check=check_output,
170+
targets=targets,
171+
compare=None,
172+
)
173+
test.run()

‎src/Model/ParticleTracking/prt.f90‎

Lines changed: 22 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,8 @@ module PrtModule
2323
use ParticleTracksModule, only: ParticleTracksType, &
2424
ParticleTrackFileType, &
2525
add_particle_event
26-
use SimModule, only: count_errors, store_error, store_error_filename
26+
use SimModule, only: count_errors, store_error, store_error_filename, &
27+
store_warning
2728
use MemoryManagerModule, only: mem_allocate
2829
use MethodModule, only: MethodType, LEVEL_FEATURE
2930
use MethodDisModule, only: MethodDisType, create_method_dis
@@ -1046,6 +1047,7 @@ subroutine prt_solve(this, isuppress_output)
10461047
use PrtPrpModule, only: PrtPrpType
10471048
use ParticleModule, only: ACTIVE, TERM_UNRELEASED, TERM_TIMEOUT
10481049
use ParticleEventModule, only: RELEASE, TERMINATE
1050+
use SimVariablesModule, only: warnmsg
10491051
! dummy
10501052
class(PrtModelType) :: this
10511053
integer(I4B), intent(in) :: isuppress_output
@@ -1099,8 +1101,25 @@ subroutine prt_solve(this, isuppress_output)
10991101
else
11001102
tmax = min(totimc + delt, particle%tstop)
11011103
end if
1102-
! Apply the tracking method until the maximum time.
1103-
call this%method%apply(particle, tmax)
1104+
! In the ordinary case tmax >= particle%ttrack always holds,
1105+
! since ttrack can't get ahead of totimc (the smaller of the
1106+
! two terms tmax is drawn from) -- so this condition failing
1107+
! really means the particle's stop time precedes its release
1108+
! time (ttrack == trelease at the point of release, which is
1109+
! the only time this branch can be reached in practice).
1110+
if (tmax >= particle%ttrack) then
1111+
! Apply the tracking method until the maximum time.
1112+
call this%method%apply(particle, tmax)
1113+
else
1114+
! STOPTIME precedes the particle's release time (see above).
1115+
write (warnmsg, '(a,g0,a,g0,a,g0,a)') &
1116+
'Particle release point ', particle%irpt, ' released at &
1117+
&time ', particle%trelease, ' but package stop time ', &
1118+
particle%tstop, ' has already passed; terminating &
1119+
&without tracking.'
1120+
call store_warning(warnmsg)
1121+
call this%method%terminate(particle, status=TERM_TIMEOUT)
1122+
end if
11041123
! If the particle timed out, terminate it.
11051124
! "Timed out" means it's still active but
11061125
! - it reached its stop time, or

0 commit comments

Comments
 (0)