Skip to content

Commit 5fe023f

Browse files
authored
feat(prt): configurable event buffering (#2789)
Previously PRT model particle events were written to output file(s) immediately. This was a deliberate choice motivated by MP7 suffering OOM errors when consolidating internal scratch files before writing the final output file. With GWF models using adaptive time stepping (ATS), or with multiple Picard iterations when PRT and GWF are in the same solution group, output would be written for discarded solve attempts. Events should only end up in output files for the final, successful solve. To gate events until a successful solve, buffer them until it's known whether they should be written. Introduce two buffering strategies, in memory and scratch file, with a new OC option SCRATCH_BUFFER enabling the latter. Memory buffering is used by default. Since the buffer could conceivably grow large for models with lots of particles or lots of action per time step (e.g., the final time step when extended tracking is on), the scratch file approach is provided to let the modeler avoid OOM errors at the price of more file IO. While this departs from informal MF6 convention to avoid scratch files in general, it's a well motivated option, not on by default. This is a pure refactor but for the new option, which I figure merits calling it a feature. It was split out of #2753. The retry misbehavior is fixed in that PR, not here.
1 parent b43e0d0 commit 5fe023f

13 files changed

Lines changed: 503 additions & 121 deletions

File tree

autotest/test_prt_voronoi1.py

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -159,7 +159,7 @@ def build_gwf_sim(name, ws, targets):
159159
}
160160

161161

162-
def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
162+
def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids, scratch_buffer=False):
163163
prt_ws = Path(prt_ws)
164164
gwf_name = get_model_name(name, "gwf")
165165
prt_name = get_model_name(name, "prt")
@@ -224,6 +224,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
224224
track_usertime=times[idx],
225225
ntracktimes=len(tracktimes) if times[idx] else None,
226226
tracktimes=[(t,) for t in tracktimes] if times[idx] else None,
227+
scratch_buffer=scratch_buffer,
227228
)
228229
gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
229230
gwf_head_file = gwf_ws / f"{gwf_name}.hds"
@@ -235,10 +236,16 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
235236
return sim
236237

237238

238-
def build_models(idx, test):
239+
def build_models(idx, test, scratch_buffer=False):
239240
gwf_sim, cell_ids = build_gwf_sim(test.name, test.workspace, test.targets)
240241
prt_sim = build_prt_sim(
241-
idx, test.name, test.workspace, test.workspace / "prt", test.targets, cell_ids
242+
idx,
243+
test.name,
244+
test.workspace,
245+
test.workspace / "prt",
246+
test.targets,
247+
cell_ids,
248+
scratch_buffer,
242249
)
243250
return gwf_sim, prt_sim
244251

@@ -357,11 +364,12 @@ def check_output(idx, test):
357364
@requires_pkg("syrupy")
358365
@pytest.mark.slow
359366
@pytest.mark.parametrize("idx, name", enumerate(cases))
360-
def test_mf6model(idx, name, function_tmpdir, targets, benchmark, plot):
367+
@pytest.mark.parametrize("scratch_buffer", [False, True], ids=["inmem", "scratch"])
368+
def test_mf6model(scratch_buffer, idx, name, function_tmpdir, targets, benchmark, plot):
361369
test = TestFramework(
362370
name=name,
363371
workspace=function_tmpdir,
364-
build=lambda t: build_models(idx, t),
372+
build=lambda t: build_models(idx, t, scratch_buffer),
365373
check=lambda t: check_output(idx, t),
366374
plot=lambda t: plot_output(idx, t) if plot else None,
367375
targets=targets,

autotest/test_prt_voronoi_perf.py

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
"""
2+
Bigger PRT benchmark on the Voronoi grid from test_prt_voronoi1.py.
3+
4+
Uses the same simple left-to-right flow field but with many more particles,
5+
generating many more track events. Meant to benchmark event buffer options,
6+
in-memory vs scratch file, at a scale where I/O overhead differences start
7+
to be measurable; very likely larger than configured here, but the runtime
8+
also grows quickly, so tune particle count and TDIS granularity as needed.
9+
"""
10+
11+
from pathlib import Path
12+
13+
import flopy
14+
import numpy as np
15+
import pytest
16+
from flopy.discretization import VertexGrid
17+
from framework import TestFramework
18+
from prt_test_utils import get_model_name
19+
from test_prt_voronoi1 import botm, build_gwf_sim, get_grid, nlay, porosity, top
20+
21+
simname = "prtvorperf"
22+
ntracktimes = 500
23+
tracktimes = list(np.linspace(0, 40000, ntracktimes))
24+
rpts = [[20, y, 0.5] for y in np.linspace(0.5, 999.5, 10000)]
25+
26+
27+
def build_prt_sim(name, gwf_ws, prt_ws, targets, scratch_buffer=False):
28+
prt_ws = Path(prt_ws)
29+
gwf_name = get_model_name(name, "gwf")
30+
prt_name = get_model_name(name, "prt")
31+
32+
grid = get_grid(prt_ws / "grid", targets)
33+
vgrid = VertexGrid(**grid.get_gridprops_vertexgrid(), nlay=1)
34+
35+
sim = flopy.mf6.MFSimulation(
36+
sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=prt_ws
37+
)
38+
flopy.mf6.ModflowTdis(sim, time_units="DAYS", perioddata=[[1.0, 1, 1.0]])
39+
prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
40+
flopy.mf6.ModflowGwfdisv(
41+
prt, nlay=nlay, **grid.get_disv_gridprops(), top=top, botm=botm
42+
)
43+
flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
44+
45+
prpdata = [
46+
(i, (0, vgrid.intersect(p[0], p[1])), p[0], p[1], p[2])
47+
for i, p in enumerate(rpts)
48+
]
49+
prp_track_file = f"{prt_name}.prp.trk"
50+
flopy.mf6.ModflowPrtprp(
51+
prt,
52+
pname="prp1",
53+
filename=f"{prt_name}_1.prp",
54+
nreleasepts=len(prpdata),
55+
packagedata=prpdata,
56+
perioddata={0: ["FIRST"]},
57+
track_filerecord=[prp_track_file],
58+
boundnames=True,
59+
stop_at_weak_sink=True,
60+
exit_solve_tolerance=1e-10,
61+
extend_tracking=True,
62+
)
63+
flopy.mf6.ModflowPrtoc(
64+
prt,
65+
pname="oc",
66+
track_exit=True,
67+
track_release=True,
68+
track_terminate=True,
69+
track_usertime=True,
70+
ntracktimes=ntracktimes,
71+
tracktimes=[(t,) for t in tracktimes],
72+
scratch_buffer=scratch_buffer,
73+
)
74+
gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
75+
gwf_head_file = gwf_ws / f"{gwf_name}.hds"
76+
flopy.mf6.ModflowPrtfmi(
77+
prt, packagedata=[("GWFHEAD", gwf_head_file), ("GWFBUDGET", gwf_budget_file)]
78+
)
79+
ems = flopy.mf6.ModflowEms(sim, pname="ems", filename=f"{prt_name}.ems")
80+
sim.register_solution_package(ems, [prt.name])
81+
return sim
82+
83+
84+
def build_models(test, scratch_buffer=False):
85+
gwf_sim, _ = build_gwf_sim(test.name, test.workspace, test.targets)
86+
prt_sim = build_prt_sim(
87+
test.name,
88+
test.workspace,
89+
test.workspace / "prt",
90+
test.targets,
91+
scratch_buffer=scratch_buffer,
92+
)
93+
return gwf_sim, prt_sim
94+
95+
96+
def check_output(test):
97+
prt_ws = test.workspace / "prt"
98+
prt_name = get_model_name(test.name, "prt")
99+
trk = prt_ws / f"{prt_name}.prp.trk"
100+
assert trk.exists() and trk.stat().st_size > 0
101+
102+
103+
@pytest.mark.skip(reason="run manually") # uncomment to run
104+
@pytest.mark.slow
105+
@pytest.mark.parametrize("scratch_buffer", [False, True], ids=["inmem", "scratch"])
106+
def test_mf6model(scratch_buffer, function_tmpdir, targets, benchmark):
107+
test = TestFramework(
108+
name=simname,
109+
workspace=function_tmpdir,
110+
build=lambda t: build_models(t, scratch_buffer),
111+
check=lambda t: check_output(t),
112+
targets=targets,
113+
compare=None,
114+
)
115+
benchmark.pedantic(test.run, rounds=1, iterations=1)

doc/ReleaseNotes/develop.toml

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,3 +88,8 @@ description = "PRT was described as configuring a single release at the beginnin
8888
section = "fixes"
8989
subsection = "basic"
9090
description = "The BMI interface's get\\_grid\\_model\\_type function worked only for numerical model types such as GWF, GWT, and GWE, but would crash for PRT models. Fix the function's implementation to work with all model types."
91+
92+
[[items]]
93+
section = "features"
94+
subsection = "model"
95+
description = "Previously PRT model particle events were written to output file(s) immediately. By default, PRT will now buffer events in memory before writing them to disk. Introduce a new option for the PRT Output Control (OC) package, SCRATCH\\_BUFFER, which enables buffering events in scratch files instead of memory. This can avoid out-of-memory errors at the price of increased file IO and runtime."

doc/mf6io/mf6ivar/dfn/prt-oc.dfn

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -292,6 +292,14 @@ mf6internal dev_dump_evtrace
292292
longname print particle tracking events
293293
description print a verbose particle tracking event trace to standard output
294294

295+
block options
296+
name scratch_buffer
297+
type keyword
298+
reader urword
299+
optional true
300+
longname buffer track events in scratch file instead of memory
301+
description keyword to stage track events in a scratch file instead of an in-memory buffer, which prevents out-of-memory errors for simulations with very many particles or events
302+
295303
# --------------------- prt oc dimensions ------------------
296304

297305
block dimensions

msvs/mf6core.vfproj

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -585,6 +585,9 @@
585585
<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ParticleEvent.f90"/>
586586
<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ParticleEvents.f90"/>
587587
<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleReleaseSchedule.f90"/>
588+
<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleTrackEventBuffer.f90"/>
589+
<File RelativePath="..\src\Solution\ParticleTracker\Particle\MemoryBuffer.f90"/>
590+
<File RelativePath="..\src\Solution\ParticleTracker\Particle\ScratchFileBuffer.f90"/>
588591
<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleTracks.f90"/>
589592
<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ReleaseEvent.f90"/>
590593
<File RelativePath="..\src\Solution\SolutionFactory.F90">

src/Idm/prt-ocidm.f90

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ module PrtOcInputModule
4040
logical :: track_timesfile = .false.
4141
logical :: timesfile = .false.
4242
logical :: dev_dump_evtrace = .false.
43+
logical :: scratch_buffer = .false.
4344
logical :: ntracktimes = .false.
4445
logical :: time = .false.
4546
logical :: saverecord = .false.
@@ -594,6 +595,25 @@ module PrtOcInputModule
594595
.false. & ! timeseries
595596
)
596597

598+
type(InputParamDefinitionType), parameter :: &
599+
prtoc_scratch_buffer = InputParamDefinitionType &
600+
( &
601+
'PRT', & ! component
602+
'OC', & ! subcomponent
603+
'OPTIONS', & ! block
604+
'SCRATCH_BUFFER', & ! tag name
605+
'SCRATCH_BUFFER', & ! fortran variable
606+
'KEYWORD', & ! type
607+
'', & ! shape
608+
'buffer track events in scratch file instead of memory', & ! longname
609+
.false., & ! required
610+
.false., & ! developmode
611+
.false., & ! multi-record
612+
.false., & ! preserve case
613+
.false., & ! layered
614+
.false. & ! timeseries
615+
)
616+
597617
type(InputParamDefinitionType), parameter :: &
598618
prtoc_ntracktimes = InputParamDefinitionType &
599619
( &
@@ -853,6 +873,7 @@ module PrtOcInputModule
853873
prtoc_track_timesfile, &
854874
prtoc_timesfile, &
855875
prtoc_dev_dump_evtrace, &
876+
prtoc_scratch_buffer, &
856877
prtoc_ntracktimes, &
857878
prtoc_time, &
858879
prtoc_saverecord, &

src/Model/ParticleTracking/prt-oc.f90

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ module PrtOcModule
3333
logical(LGP), pointer :: trackdropped => null() !< whether to track drops to water table
3434
integer(I4B), pointer :: ntracktimes => null() !< number of user-specified tracking times
3535
logical(LGP), pointer :: dump_event_trace => null() !< whether to dump event trace for debugging
36+
logical(LGP), pointer :: scratch_buffer => null() !< whether to use scratch file instead of memory for event buffering
3637
type(TimeSelectType), pointer :: tracktimes !< user-specified tracking times
3738

3839
contains
@@ -78,6 +79,7 @@ subroutine prt_oc_allocate_scalars(this, name_model, input_mempath)
7879
allocate (this%name_model)
7980
allocate (this%input_fname)
8081
call mem_allocate(this%dump_event_trace, 'DUMP_EVENT_TRACE', this%memoryPath)
82+
call mem_allocate(this%scratch_buffer, 'SCRATCH_BUFFER', this%memoryPath)
8183
call mem_allocate(this%inunit, 'INUNIT', this%memoryPath)
8284
call mem_allocate(this%iout, 'IOUT', this%memoryPath)
8385
call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
@@ -101,6 +103,7 @@ subroutine prt_oc_allocate_scalars(this, name_model, input_mempath)
101103
this%input_mempath = input_mempath
102104
this%input_fname = ''
103105
this%dump_event_trace = .false.
106+
this%scratch_buffer = .false.
104107
this%inunit = 0
105108
this%iout = 0
106109
this%ibudcsv = 0
@@ -178,6 +181,7 @@ subroutine prt_oc_da(this)
178181

179182
deallocate (this%name_model)
180183
call mem_deallocate(this%dump_event_trace)
184+
call mem_deallocate(this%scratch_buffer)
181185
call mem_deallocate(this%inunit)
182186
call mem_deallocate(this%iout)
183187
call mem_deallocate(this%ibudcsv)
@@ -251,6 +255,8 @@ subroutine prt_oc_source_options(this)
251255
found%track_usertime)
252256
call mem_set_value(evinput, 'DEV_DUMP_EVTRACE', this%input_mempath, &
253257
found%dev_dump_evtrace)
258+
call mem_set_value(evinput, 'SCRATCH_BUFFER', this%input_mempath, &
259+
found%scratch_buffer)
254260

255261
if (found%track_release) this%trackrelease = .true.
256262
if (found%track_exit) this%trackfeatexit = .true.
@@ -261,6 +267,13 @@ subroutine prt_oc_source_options(this)
261267
if (found%track_weaksink) this%trackweaksink = .true.
262268
if (found%track_usertime) this%trackusertime = .true.
263269
if (found%dev_dump_evtrace) this%dump_event_trace = .true.
270+
if (found%scratch_buffer) this%scratch_buffer = .true.
271+
272+
if (this%scratch_buffer) then
273+
write (this%iout, '(4x,a)') 'TRACK EVENT BUFFER: SCRATCH FILE'
274+
else
275+
write (this%iout, '(4x,a)') 'TRACK EVENT BUFFER: MEMORY'
276+
end if
264277

265278
! default to all events
266279
if (.not. (found%track_release .or. &

src/Model/ParticleTracking/prt.f90

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ module PrtModule
2222
use ParticleEventsModule, only: ParticleEventDispatcherType, handle_event
2323
use ParticleTracksModule, only: ParticleTracksType, &
2424
ParticleTrackFileType, &
25-
write_particle_event
25+
add_particle_event
2626
use SimModule, only: count_errors, store_error, store_error_filename
2727
use MemoryManagerModule, only: mem_allocate
2828
use MethodModule, only: MethodType, LEVEL_FEATURE
@@ -53,7 +53,7 @@ module PrtModule
5353
type(MethodDisType), pointer :: method_dis => null() ! DIS tracking method
5454
type(MethodDisvType), pointer :: method_disv => null() ! DISV tracking method
5555
type(ParticleEventDispatcherType), pointer :: events => null() ! event dispatcher
56-
class(ParticleTracksType), pointer :: tracks ! track output manager
56+
type(ParticleTracksType), pointer :: tracks ! track output manager
5757
integer(I4B), pointer :: infmi => null() ! unit number FMI
5858
integer(I4B), pointer :: inmip => null() ! unit number MIP
5959
integer(I4B), pointer :: inmvt => null() ! unit number MVT
@@ -262,6 +262,9 @@ subroutine prt_ar(this)
262262
call this%oc%oc_ar(this%dis, DHNOFLO)
263263
call this%budget%set_ibudcsv(this%oc%ibudcsv)
264264

265+
! Initialize the event buffer (memory or scratch file per OC option)
266+
call this%tracks%init_buffer(this%oc%scratch_buffer)
267+
265268
! Select tracking events
266269
call this%tracks%select_events( &
267270
this%oc%trackrelease, &
@@ -331,7 +334,7 @@ subroutine prt_ar(this)
331334

332335
! Subscribe particle track output manager to events
333336
p => this%tracks
334-
call this%events%subscribe(write_particle_event, p)
337+
call this%events%subscribe(add_particle_event, p)
335338

336339
! Set verbose tracing if requested
337340
if (this%oc%dump_event_trace) this%tracks%iout = 0
@@ -573,7 +576,8 @@ subroutine prt_ot(this)
573576
integer(I4B) :: ibudfl
574577
integer(I4B) :: ipflag
575578

576-
! Note: particle tracking output is handled elsewhere
579+
! Flush buffered events to disk
580+
call this%tracks%flush_buffer()
577581

578582
! Set write and print flags
579583
idvsave = 0

0 commit comments

Comments
 (0)