feat(prt): configurable event buffering

autotest/test_prt_voronoi1.py

@@ -159,7 +159,7 @@ def build_gwf_sim(name, ws, targets):
     }
 
 
-def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
+def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids, scratch_buffer=False):
     prt_ws = Path(prt_ws)
     gwf_name = get_model_name(name, "gwf")
     prt_name = get_model_name(name, "prt")
@@ -224,6 +224,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
         track_usertime=times[idx],
         ntracktimes=len(tracktimes) if times[idx] else None,
         tracktimes=[(t,) for t in tracktimes] if times[idx] else None,
+        scratch_buffer=scratch_buffer,
     )
     gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
     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):
     return sim
 
 
-def build_models(idx, test):
+def build_models(idx, test, scratch_buffer=False):
     gwf_sim, cell_ids = build_gwf_sim(test.name, test.workspace, test.targets)
     prt_sim = build_prt_sim(
-        idx, test.name, test.workspace, test.workspace / "prt", test.targets, cell_ids
+        idx,
+        test.name,
+        test.workspace,
+        test.workspace / "prt",
+        test.targets,
+        cell_ids,
+        scratch_buffer,
     )
     return gwf_sim, prt_sim
 
@@ -357,11 +364,12 @@ def check_output(idx, test):
 @requires_pkg("syrupy")
 @pytest.mark.slow
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets, benchmark, plot):
+@pytest.mark.parametrize("scratch_buffer", [False, True], ids=["inmem", "scratch"])
+def test_mf6model(scratch_buffer, idx, name, function_tmpdir, targets, benchmark, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
-        build=lambda t: build_models(idx, t),
+        build=lambda t: build_models(idx, t, scratch_buffer),
         check=lambda t: check_output(idx, t),
         plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,

autotest/test_prt_voronoi_perf.py

@@ -0,0 +1,115 @@
+"""
+Bigger PRT benchmark on the Voronoi grid from test_prt_voronoi1.py.
+
+Uses the same simple left-to-right flow field but with many more particles,
+generating many more track events. Meant to benchmark event buffer options,
+in-memory vs scratch file, at a scale where I/O overhead differences start
+to be measurable; very likely larger than configured here, but the runtime
+also grows quickly, so tune particle count and TDIS granularity as needed.
+"""
+
+from pathlib import Path
+
+import flopy
+import numpy as np
+import pytest
+from flopy.discretization import VertexGrid
+from framework import TestFramework
+from prt_test_utils import get_model_name
+from test_prt_voronoi1 import botm, build_gwf_sim, get_grid, nlay, porosity, top
+
+simname = "prtvorperf"
+ntracktimes = 500
+tracktimes = list(np.linspace(0, 40000, ntracktimes))
+rpts = [[20, y, 0.5] for y in np.linspace(0.5, 999.5, 10000)]
+
+
+def build_prt_sim(name, gwf_ws, prt_ws, targets, scratch_buffer=False):
+    prt_ws = Path(prt_ws)
+    gwf_name = get_model_name(name, "gwf")
+    prt_name = get_model_name(name, "prt")
+
+    grid = get_grid(prt_ws / "grid", targets)
+    vgrid = VertexGrid(**grid.get_gridprops_vertexgrid(), nlay=1)
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=prt_ws
+    )
+    flopy.mf6.ModflowTdis(sim, time_units="DAYS", perioddata=[[1.0, 1, 1.0]])
+    prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
+    flopy.mf6.ModflowGwfdisv(
+        prt, nlay=nlay, **grid.get_disv_gridprops(), top=top, botm=botm
+    )
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+
+    prpdata = [
+        (i, (0, vgrid.intersect(p[0], p[1])), p[0], p[1], p[2])
+        for i, p in enumerate(rpts)
+    ]
+    prp_track_file = f"{prt_name}.prp.trk"
+    flopy.mf6.ModflowPrtprp(
+        prt,
+        pname="prp1",
+        filename=f"{prt_name}_1.prp",
+        nreleasepts=len(prpdata),
+        packagedata=prpdata,
+        perioddata={0: ["FIRST"]},
+        track_filerecord=[prp_track_file],
+        boundnames=True,
+        stop_at_weak_sink=True,
+        exit_solve_tolerance=1e-10,
+        extend_tracking=True,
+    )
+    flopy.mf6.ModflowPrtoc(
+        prt,
+        pname="oc",
+        track_exit=True,
+        track_release=True,
+        track_terminate=True,
+        track_usertime=True,
+        ntracktimes=ntracktimes,
+        tracktimes=[(t,) for t in tracktimes],
+        scratch_buffer=scratch_buffer,
+    )
+    gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
+    gwf_head_file = gwf_ws / f"{gwf_name}.hds"
+    flopy.mf6.ModflowPrtfmi(
+        prt, packagedata=[("GWFHEAD", gwf_head_file), ("GWFBUDGET", gwf_budget_file)]
+    )
+    ems = flopy.mf6.ModflowEms(sim, pname="ems", filename=f"{prt_name}.ems")
+    sim.register_solution_package(ems, [prt.name])
+    return sim
+
+
+def build_models(test, scratch_buffer=False):
+    gwf_sim, _ = build_gwf_sim(test.name, test.workspace, test.targets)
+    prt_sim = build_prt_sim(
+        test.name,
+        test.workspace,
+        test.workspace / "prt",
+        test.targets,
+        scratch_buffer=scratch_buffer,
+    )
+    return gwf_sim, prt_sim
+
+
+def check_output(test):
+    prt_ws = test.workspace / "prt"
+    prt_name = get_model_name(test.name, "prt")
+    trk = prt_ws / f"{prt_name}.prp.trk"
+    assert trk.exists() and trk.stat().st_size > 0
+
+
+@pytest.mark.skip(reason="run manually")  # uncomment to run
+@pytest.mark.slow
+@pytest.mark.parametrize("scratch_buffer", [False, True], ids=["inmem", "scratch"])
+def test_mf6model(scratch_buffer, function_tmpdir, targets, benchmark):
+    test = TestFramework(
+        name=simname,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(t, scratch_buffer),
+        check=lambda t: check_output(t),
+        targets=targets,
+        compare=None,
+    )
+    benchmark.pedantic(test.run, rounds=1, iterations=1)

doc/ReleaseNotes/develop.toml

@@ -88,3 +88,8 @@ description = "PRT was described as configuring a single release at the beginnin
 section = "fixes"
 subsection = "basic"
 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."
+
+[[items]]
+section = "features"
+subsection = "model"
+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

@@ -292,6 +292,14 @@ mf6internal dev_dump_evtrace
 longname print particle tracking events
 description print a verbose particle tracking event trace to standard output
 
+block options
+name scratch_buffer
+type keyword
+reader urword
+optional true
+longname buffer track events in scratch file instead of memory
+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
+
 # --------------------- prt oc dimensions ------------------
 
 block dimensions

msvs/mf6core.vfproj

@@ -585,6 +585,9 @@
 				<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ParticleEvent.f90"/>
 				<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ParticleEvents.f90"/>
 				<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleReleaseSchedule.f90"/>
+				<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleTrackEventBuffer.f90"/>
+				<File RelativePath="..\src\Solution\ParticleTracker\Particle\MemoryBuffer.f90"/>
+				<File RelativePath="..\src\Solution\ParticleTracker\Particle\ScratchFileBuffer.f90"/>
 				<File RelativePath="..\src\Solution\ParticleTracker\Particle\ParticleTracks.f90"/>
 				<File RelativePath="..\src\Solution\ParticleTracker\Particle\Events\ReleaseEvent.f90"/>
 				<File RelativePath="..\src\Solution\SolutionFactory.F90">

src/Idm/prt-ocidm.f90

@@ -40,6 +40,7 @@ module PrtOcInputModule
     logical :: track_timesfile = .false.
     logical :: timesfile = .false.
     logical :: dev_dump_evtrace = .false.
+    logical :: scratch_buffer = .false.
     logical :: ntracktimes = .false.
     logical :: time = .false.
     logical :: saverecord = .false.
@@ -594,6 +595,25 @@ module PrtOcInputModule
     .false. & ! timeseries
     )
 
+  type(InputParamDefinitionType), parameter :: &
+    prtoc_scratch_buffer = InputParamDefinitionType &
+    ( &
+    'PRT', & ! component
+    'OC', & ! subcomponent
+    'OPTIONS', & ! block
+    'SCRATCH_BUFFER', & ! tag name
+    'SCRATCH_BUFFER', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    'buffer track events in scratch file instead of memory', & ! longname
+    .false., & ! required
+    .false., & ! developmode
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
   type(InputParamDefinitionType), parameter :: &
     prtoc_ntracktimes = InputParamDefinitionType &
     ( &
@@ -853,6 +873,7 @@ module PrtOcInputModule
     prtoc_track_timesfile, &
     prtoc_timesfile, &
     prtoc_dev_dump_evtrace, &
+    prtoc_scratch_buffer, &
     prtoc_ntracktimes, &
     prtoc_time, &
     prtoc_saverecord, &

src/Model/ParticleTracking/prt-oc.f90

@@ -33,6 +33,7 @@ module PrtOcModule
     logical(LGP), pointer :: trackdropped => null() !< whether to track drops to water table
     integer(I4B), pointer :: ntracktimes => null() !< number of user-specified tracking times
     logical(LGP), pointer :: dump_event_trace => null() !< whether to dump event trace for debugging
+    logical(LGP), pointer :: scratch_buffer => null() !< whether to use scratch file instead of memory for event buffering
     type(TimeSelectType), pointer :: tracktimes !< user-specified tracking times
 
   contains
@@ -78,6 +79,7 @@ subroutine prt_oc_allocate_scalars(this, name_model, input_mempath)
     allocate (this%name_model)
     allocate (this%input_fname)
     call mem_allocate(this%dump_event_trace, 'DUMP_EVENT_TRACE', this%memoryPath)
+    call mem_allocate(this%scratch_buffer, 'SCRATCH_BUFFER', this%memoryPath)
     call mem_allocate(this%inunit, 'INUNIT', this%memoryPath)
     call mem_allocate(this%iout, 'IOUT', this%memoryPath)
     call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
@@ -101,6 +103,7 @@ subroutine prt_oc_allocate_scalars(this, name_model, input_mempath)
     this%input_mempath = input_mempath
     this%input_fname = ''
     this%dump_event_trace = .false.
+    this%scratch_buffer = .false.
     this%inunit = 0
     this%iout = 0
     this%ibudcsv = 0
@@ -178,6 +181,7 @@ subroutine prt_oc_da(this)
 
     deallocate (this%name_model)
     call mem_deallocate(this%dump_event_trace)
+    call mem_deallocate(this%scratch_buffer)
     call mem_deallocate(this%inunit)
     call mem_deallocate(this%iout)
     call mem_deallocate(this%ibudcsv)
@@ -251,6 +255,8 @@ subroutine prt_oc_source_options(this)
                        found%track_usertime)
     call mem_set_value(evinput, 'DEV_DUMP_EVTRACE', this%input_mempath, &
                        found%dev_dump_evtrace)
+    call mem_set_value(evinput, 'SCRATCH_BUFFER', this%input_mempath, &
+                       found%scratch_buffer)
 
     if (found%track_release) this%trackrelease = .true.
     if (found%track_exit) this%trackfeatexit = .true.
@@ -261,6 +267,13 @@ subroutine prt_oc_source_options(this)
     if (found%track_weaksink) this%trackweaksink = .true.
     if (found%track_usertime) this%trackusertime = .true.
     if (found%dev_dump_evtrace) this%dump_event_trace = .true.
+    if (found%scratch_buffer) this%scratch_buffer = .true.
+
+    if (this%scratch_buffer) then
+      write (this%iout, '(4x,a)') 'TRACK EVENT BUFFER: SCRATCH FILE'
+    else
+      write (this%iout, '(4x,a)') 'TRACK EVENT BUFFER: MEMORY'
+    end if
 
     ! default to all events
     if (.not. (found%track_release .or. &

src/Model/ParticleTracking/prt.f90

@@ -22,7 +22,7 @@ module PrtModule
   use ParticleEventsModule, only: ParticleEventDispatcherType, handle_event
   use ParticleTracksModule, only: ParticleTracksType, &
                                   ParticleTrackFileType, &
-                                  write_particle_event
+                                  add_particle_event
   use SimModule, only: count_errors, store_error, store_error_filename
   use MemoryManagerModule, only: mem_allocate
   use MethodModule, only: MethodType, LEVEL_FEATURE
@@ -53,7 +53,7 @@ module PrtModule
     type(MethodDisType), pointer :: method_dis => null() ! DIS tracking method
     type(MethodDisvType), pointer :: method_disv => null() ! DISV tracking method
     type(ParticleEventDispatcherType), pointer :: events => null() ! event dispatcher
-    class(ParticleTracksType), pointer :: tracks ! track output manager
+    type(ParticleTracksType), pointer :: tracks ! track output manager
     integer(I4B), pointer :: infmi => null() ! unit number FMI
     integer(I4B), pointer :: inmip => null() ! unit number MIP
     integer(I4B), pointer :: inmvt => null() ! unit number MVT
@@ -262,6 +262,9 @@ subroutine prt_ar(this)
     call this%oc%oc_ar(this%dis, DHNOFLO)
     call this%budget%set_ibudcsv(this%oc%ibudcsv)
 
+    ! Initialize the event buffer (memory or scratch file per OC option)
+    call this%tracks%init_buffer(this%oc%scratch_buffer)
+
     ! Select tracking events
     call this%tracks%select_events( &
       this%oc%trackrelease, &
@@ -331,7 +334,7 @@ subroutine prt_ar(this)
 
     ! Subscribe particle track output manager to events
     p => this%tracks
-    call this%events%subscribe(write_particle_event, p)
+    call this%events%subscribe(add_particle_event, p)
 
     ! Set verbose tracing if requested
     if (this%oc%dump_event_trace) this%tracks%iout = 0
@@ -573,7 +576,8 @@ subroutine prt_ot(this)
     integer(I4B) :: ibudfl
     integer(I4B) :: ipflag
 
-    ! Note: particle tracking output is handled elsewhere
+    ! Flush buffered events to disk
+    call this%tracks%flush_buffer()
 
     ! Set write and print flags
     idvsave = 0