diff --git a/autotest/test_prt_voronoi1.py b/autotest/test_prt_voronoi1.py index a36eee6c4b4..f4d9a4549b5 100644 --- a/autotest/test_prt_voronoi1.py +++ b/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, diff --git a/autotest/test_prt_voronoi_perf.py b/autotest/test_prt_voronoi_perf.py new file mode 100644 index 00000000000..28632a59b48 --- /dev/null +++ b/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) diff --git a/doc/ReleaseNotes/develop.toml b/doc/ReleaseNotes/develop.toml index 7e316e53e65..0b5a6b72014 100644 --- a/doc/ReleaseNotes/develop.toml +++ b/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." diff --git a/doc/mf6io/mf6ivar/dfn/prt-oc.dfn b/doc/mf6io/mf6ivar/dfn/prt-oc.dfn index 4bb51c2ae41..03a170d09ad 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-oc.dfn +++ b/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 diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 742f5bbfbba..53a3d5e1063 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -585,6 +585,9 @@ + + + diff --git a/src/Idm/prt-ocidm.f90 b/src/Idm/prt-ocidm.f90 index c375832c3c3..b72479da887 100644 --- a/src/Idm/prt-ocidm.f90 +++ b/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, & diff --git a/src/Model/ParticleTracking/prt-oc.f90 b/src/Model/ParticleTracking/prt-oc.f90 index fd43446f3dd..09eeebacc2d 100644 --- a/src/Model/ParticleTracking/prt-oc.f90 +++ b/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. & diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index 3d19d00b2c3..73eb1ac6549 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/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 diff --git a/src/Solution/ParticleTracker/Particle/MemoryBuffer.f90 b/src/Solution/ParticleTracker/Particle/MemoryBuffer.f90 new file mode 100644 index 00000000000..661e0282e9b --- /dev/null +++ b/src/Solution/ParticleTracker/Particle/MemoryBuffer.f90 @@ -0,0 +1,71 @@ +!> @brief In-memory particle event buffer module. +!< +module MemoryBufferModule + + use KindModule, only: I4B + use ParticleTrackEventBufferModule, only: ParticleTrackEventBufferType, & + ParticleTrackRecordType, & + ParticleTrackFileType, & + save_record + + implicit none + private + public :: MemoryBufferType + + !> @brief In-memory particle event buffer. Records are held in a + !! dynamically growing array that doubles in capacity as needed. + type, extends(ParticleTrackEventBufferType) :: MemoryBufferType + type(ParticleTrackRecordType), allocatable :: records(:) !< buffer + contains + procedure :: append => memory_append + procedure :: flush => memory_flush + procedure :: discard => memory_discard + procedure :: destroy => memory_destroy + end type MemoryBufferType + +contains + + subroutine memory_append(this, rec) + class(MemoryBufferType) :: this + type(ParticleTrackRecordType), intent(in) :: rec + type(ParticleTrackRecordType), allocatable :: tmp(:) + + if (.not. allocated(this%records)) then + allocate (this%records(64)) + else if (this%nrecords == size(this%records)) then + allocate (tmp(size(this%records) * 2)) + tmp(1:this%nrecords) = this%records + call move_alloc(tmp, this%records) + end if + this%nrecords = this%nrecords + 1 + this%records(this%nrecords) = rec + end subroutine memory_append + + subroutine memory_flush(this, files) + class(MemoryBufferType) :: this + type(ParticleTrackFileType), intent(in) :: files(:) + integer(I4B) :: n, i + type(ParticleTrackRecordType) :: rec + + do n = 1, this%nrecords + rec = this%records(n) + do i = 1, size(files) + if (files(i)%iun > 0 .and. & + (files(i)%iprp == -1 .or. files(i)%iprp == rec%iprp)) & + call save_record(files(i)%iun, rec, files(i)%csv) + end do + end do + this%nrecords = 0 + end subroutine memory_flush + + subroutine memory_discard(this) + class(MemoryBufferType) :: this + this%nrecords = 0 + end subroutine memory_discard + + subroutine memory_destroy(this) + class(MemoryBufferType) :: this + if (allocated(this%records)) deallocate (this%records) + end subroutine memory_destroy + +end module MemoryBufferModule diff --git a/src/Solution/ParticleTracker/Particle/ParticleTrackEventBuffer.f90 b/src/Solution/ParticleTracker/Particle/ParticleTrackEventBuffer.f90 new file mode 100644 index 00000000000..1de58a45395 --- /dev/null +++ b/src/Solution/ParticleTracker/Particle/ParticleTrackEventBuffer.f90 @@ -0,0 +1,90 @@ +!> @brief Particle track event buffer module. +!! +!! Defines the abstract ParticleTrackEventBufferType, the ParticleTrackRecordType +!! and ParticleTrackFileType data types, and the save_record helper used +!! by concrete buffer implementations. +!< +module ParticleTrackEventBufferModule + + use KindModule, only: DP, I4B, LGP + use ErrorUtilModule, only: pstop + + implicit none + + !> @brief Flat record of a particle track event. + type :: ParticleTrackRecordType + integer(I4B) :: kper, kstp, imdl, iprp, irpt, ilay, icu, izone + integer(I4B) :: istatus, ireason + real(DP) :: trelease, ttrack, x, y, z + character(len=40) :: name + end type ParticleTrackRecordType + + !> @brief Output file containing all or some particle pathlines. + !! + !! Can be associated with a particle release point (PRP) package + !! or with an entire model, and can be binary or text (CSV). + !< + type :: ParticleTrackFileType + private + integer(I4B), public :: iun = 0 !< file unit number + logical(LGP), public :: csv = .false. !< whether the file is binary or CSV + integer(I4B), public :: iprp = -1 !< -1 is model-level file, 0 is exchange PRP + end type ParticleTrackFileType + + !> @brief Event buffering strategy + type, abstract :: ParticleTrackEventBufferType + integer(I4B) :: nrecords = 0 !< number of records stored + contains + procedure :: init => buffer_init !< open/allocate resources + procedure(buffer_append), deferred :: append !< buffer one record + procedure(buffer_flush), deferred :: flush !< write buffered records, reset + procedure(buffer_simple), deferred :: discard !< reset without writing + procedure(buffer_simple), deferred :: destroy !< release resources + end type ParticleTrackEventBufferType + + abstract interface + subroutine buffer_append(this, rec) + import ParticleTrackEventBufferType, ParticleTrackRecordType + class(ParticleTrackEventBufferType) :: this + type(ParticleTrackRecordType), intent(in) :: rec + end subroutine buffer_append + + subroutine buffer_flush(this, files) + import ParticleTrackEventBufferType, ParticleTrackFileType + class(ParticleTrackEventBufferType) :: this + type(ParticleTrackFileType), intent(in) :: files(:) + end subroutine buffer_flush + + ! Shared interface for discard and destroy: both take only `this`. + subroutine buffer_simple(this) + import ParticleTrackEventBufferType + class(ParticleTrackEventBufferType) :: this + end subroutine buffer_simple + end interface + +contains + + subroutine buffer_init(this) + class(ParticleTrackEventBufferType) :: this + end subroutine buffer_init + + !> @brief Save an event record to a binary or CSV file. + subroutine save_record(iun, rec, csv) + integer(I4B), intent(in) :: iun + type(ParticleTrackRecordType), intent(in) :: rec + logical(LGP), intent(in) :: csv + + if (csv) then + write (iun, '(*(G0,:,","))') & + rec%kper, rec%kstp, rec%imdl, rec%iprp, rec%irpt, & + rec%ilay, rec%icu, rec%izone, rec%istatus, rec%ireason, & + rec%trelease, rec%ttrack, rec%x, rec%y, rec%z, trim(rec%name) + else + write (iun) & + rec%kper, rec%kstp, rec%imdl, rec%iprp, rec%irpt, & + rec%ilay, rec%icu, rec%izone, rec%istatus, rec%ireason, & + rec%trelease, rec%ttrack, rec%x, rec%y, rec%z, rec%name + end if + end subroutine save_record + +end module ParticleTrackEventBufferModule diff --git a/src/Solution/ParticleTracker/Particle/ParticleTracks.f90 b/src/Solution/ParticleTracker/Particle/ParticleTracks.f90 index 03b2d228c99..018ee9c8103 100644 --- a/src/Solution/ParticleTracker/Particle/ParticleTracks.f90 +++ b/src/Solution/ParticleTracker/Particle/ParticleTracks.f90 @@ -13,12 +13,14 @@ !! them to one or more track files, binary or CSV, and for logging the !! events if requested. Each track file is associated with either a PRP !! package or with the full PRT model (there may only be 1 such latter). +!! +!! Events can be buffered in memory or in a scratch file, and in either +!! case are flushed to disk only when a time step successfully finishes. !< module ParticleTracksModule use KindModule, only: DP, I4B, LGP use ErrorUtilModule, only: pstop - use ConstantsModule, only: DZERO, DONE, DPIO180 use ParticleModule, only: ParticleType, ACTIVE use ParticleEventModule, only: ParticleEventType use ReleaseEventModule, only: ReleaseEventType @@ -32,13 +34,14 @@ module ParticleTracksModule use ParticleEventsModule, only: ParticleEventDispatcherType use BaseDisModule, only: DisBaseType use GeomUtilModule, only: transform + use ParticleTrackEventBufferModule + use MemoryBufferModule, only: MemoryBufferType + use ScratchFileBufferModule, only: ScratchFileBufferType implicit none - public :: ParticleTrackFileType, & - ParticleTracksType, & + public :: ParticleTracksType, & ParticleTrackEventSelectionType, & - write_particle_event - private :: save_event + add_particle_event character(len=*), parameter, public :: TRACKHEADER = & 'kper,kstp,imdl,iprp,irpt,ilay,icell,izone,& @@ -48,18 +51,6 @@ module ParticleTracksModule ' @brief Output file containing all or some particle pathlines. - !! - !! Can be associated with a particle release point (PRP) package - !! or with an entire model, and can be binary or text (CSV). - !< - type :: ParticleTrackFileType - private - integer(I4B), public :: iun = 0 !< file unit number - logical(LGP), public :: csv = .false. !< whether the file is binary or CSV - integer(I4B), public :: iprp = -1 !< -1 is model-level file, 0 is exchange PRP - end type ParticleTrackFileType - !> @brief Selection of particle events. type :: ParticleTrackEventSelectionType logical(LGP) :: release !< track release events @@ -76,41 +67,41 @@ module ParticleTracksModule !! to files. One output unit can be configured for printing. Multiple files !! can be configured for writing, with each file optionally associated with !! a PRP package or with the full model. Events can be filtered by type, so - !! that only certain event types are printed or written to files. + !! that only certain event types are printed or written to files. Particle + !! events are buffered in memory or in a scratch file, flushed to disk only + !! when the time step is successfully solved for the last time (there may + !! be multiple solves per time step, depending on ATS and Picard options). !< type :: ParticleTracksType private - integer(I4B), public :: iout = -1 !< log file unit + integer(I4B), public :: iout = -1 !< log file unit integer(I4B), public :: ntrackfiles !< number of track files type(ParticleTrackFileType), public, allocatable :: files(:) !< track files type(ParticleTrackEventSelectionType), public :: selected !< event selection + class(ParticleTrackEventBufferType), allocatable :: buffer !< event buffer contains procedure, public :: init_file + procedure, public :: init_buffer procedure, public :: is_selected procedure, public :: select_events + procedure, public :: buffer_event + procedure, public :: flush_buffer + procedure, public :: discard_buffer procedure, public :: destroy procedure :: expand_files - procedure :: should_save - procedure :: should_print end type ParticleTracksType contains - !> @brief Initialize a binary or CSV file. + !> @brief Initialize a binary or CSV output file. subroutine init_file(this, iun, csv, iprp) - ! dummy class(ParticleTracksType) :: this integer(I4B), intent(in) :: iun logical(LGP), intent(in), optional :: csv integer(I4B), intent(in), optional :: iprp - ! local type(ParticleTrackFileType), pointer :: file - if (.not. allocated(this%files)) then - allocate (this%files(1)) - else - call this%expand_files(increment=1) - end if + call this%expand_files(increment=1) allocate (file) file%iun = iun @@ -120,21 +111,30 @@ subroutine init_file(this, iun, csv, iprp) this%files(this%ntrackfiles) = file end subroutine init_file + !> @brief Initialize the event buffer strategy + subroutine init_buffer(this, scratch) + class(ParticleTracksType) :: this + logical(LGP), intent(in) :: scratch + + if (scratch) then + allocate (ScratchFileBufferType :: this%buffer) + else + allocate (MemoryBufferType :: this%buffer) + end if + call this%buffer%init() + end subroutine init_buffer + !> @brief Destroy the particle track manager. subroutine destroy(this) class(ParticleTracksType) :: this - if (allocated(this%files)) deallocate (this%files) + call this%buffer%destroy() end subroutine destroy !> @brief Grow the array of track files. subroutine expand_files(this, increment) - ! dummy class(ParticleTracksType) :: this integer(I4B), optional, intent(in) :: increment - ! local - integer(I4B) :: inclocal - integer(I4B) :: isize - integer(I4B) :: newsize + integer(I4B) :: inclocal, isize, newsize type(ParticleTrackFileType), allocatable, dimension(:) :: temp if (present(increment)) then @@ -210,99 +210,59 @@ logical function is_selected(this, event) result(selected) call pstop(1, "unknown event type") selected = .false. end select - end function is_selected - !> @brief Check whether a particle belongs in a given file i.e. - !! if the file is enabled and its group matches the particle's. - logical function should_save(this, particle, file) result(save) - class(ParticleTracksType), intent(inout) :: this - type(ParticleType), pointer, intent(in) :: particle - type(ParticleTrackFileType), intent(in) :: file - save = (file%iun > 0 .and. & - (file%iprp == -1 .or. file%iprp == particle%iprp)) - end function should_save - - !> @brief Save an event to a binary or CSV file. - subroutine save_event(iun, particle, event, csv) - ! dummy - integer(I4B), intent(in) :: iun + !> @brief Buffer an event for deferred write. + subroutine buffer_event(this, particle, event) + class(ParticleTracksType) :: this type(ParticleType), pointer, intent(in) :: particle class(ParticleEventType), pointer, intent(in) :: event - logical(LGP), intent(in) :: csv + type(ParticleTrackRecordType) :: rec - if (csv) then - write (iun, '(*(G0,:,","))') & - event%kper, & - event%kstp, & - event%imdl, & - event%iprp, & - event%irpt, & - event%ilay, & - event%icu, & - event%izone, & - event%istatus, & - event%get_code(), & - event%trelease, & - event%ttrack, & - event%x, & - event%y, & - event%z, & - trim(adjustl(particle%name)) - else - write (iun) & - event%kper, & - event%kstp, & - event%imdl, & - event%iprp, & - event%irpt, & - event%ilay, & - event%icu, & - event%izone, & - event%istatus, & - event%get_code(), & - event%trelease, & - event%ttrack, & - event%x, & - event%y, & - event%z, & - particle%name - end if - end subroutine save_event + rec%kper = event%kper; rec%kstp = event%kstp + rec%imdl = event%imdl; rec%iprp = event%iprp + rec%irpt = event%irpt; rec%ilay = event%ilay + rec%icu = event%icu; rec%izone = event%izone + rec%istatus = event%istatus; rec%ireason = event%get_code() + rec%trelease = event%trelease; rec%ttrack = event%ttrack + rec%x = event%x; rec%y = event%y; rec%z = event%z + rec%name = trim(adjustl(particle%name)) - !> @brief Is the output unit valid? - logical function should_print(this) - class(ParticleTracksType), intent(inout) :: this - should_print = this%iout >= 0 - end function should_print + call this%buffer%append(rec) + end subroutine buffer_event + + !> @brief Flush the event buffer to disk. + subroutine flush_buffer(this) + class(ParticleTracksType) :: this + call this%buffer%flush(this%files) + end subroutine flush_buffer + + !> @brief Discard buffered events without writing. + subroutine discard_buffer(this) + class(ParticleTracksType) :: this + call this%buffer%discard() + end subroutine discard_buffer - !> @brief Write a particle event to files for which the particle - !! is eligible, and print the event to output unit if requested. - !! This function is the module's main entry point. It should be - !! subscribed as an event handler to particle event dispatchers. - function write_particle_event(context, particle, event) result(handled) - ! dummy + !> @brief Add a particle event to be written to eligible + !! files and printed to an output file unit if requested. + !! This function should be subscribed as an event handler + !! to particle event dispatchers. Events are buffered, to + !! be written to output files upon successful completion + !! of a time step when the framework OT hook is executed. + function add_particle_event(context, particle, event) result(handled) class(*), pointer :: context type(ParticleType), pointer, intent(inout) :: particle class(ParticleEventType), pointer, intent(in) :: event logical(LGP) :: handled - ! local - integer(I4B) :: i - type(ParticleTrackFileType) :: file select type (context) type is (ParticleTracksType) - if (context%should_print()) & + if (context%iout >= 0) & call event%log(context%iout) - if (context%is_selected(event)) then - do i = 1, context%ntrackfiles - file = context%files(i) - if (context%should_save(particle, file)) & - call save_event(file%iun, particle, event, csv=file%csv) - end do - end if + if (context%is_selected(event)) & + call context%buffer_event(particle, event) handled = .true. end select - end function write_particle_event + end function add_particle_event end module ParticleTracksModule diff --git a/src/Solution/ParticleTracker/Particle/ScratchFileBuffer.f90 b/src/Solution/ParticleTracker/Particle/ScratchFileBuffer.f90 new file mode 100644 index 00000000000..f233db70d2d --- /dev/null +++ b/src/Solution/ParticleTracker/Particle/ScratchFileBuffer.f90 @@ -0,0 +1,81 @@ +!> @brief Scratch-file particle event buffer module. +!< +module ScratchFileBufferModule + + use KindModule, only: I4B + use ErrorUtilModule, only: pstop + use ParticleTrackEventBufferModule, only: ParticleTrackEventBufferType, & + ParticleTrackRecordType, & + ParticleTrackFileType, & + save_record + + implicit none + private + public :: ScratchFileBufferType + + !> @brief Scratch-file particle event buffer. Records are written to + !! an unformatted sequential scratch file that is rewound on discard + !! and read back sequentially on flush. + type, extends(ParticleTrackEventBufferType) :: ScratchFileBufferType + integer(I4B) :: iun = 0 !< scratch file unit + contains + procedure :: init => scratch_init + procedure :: append => scratch_append + procedure :: flush => scratch_flush + procedure :: discard => scratch_discard + procedure :: destroy => scratch_destroy + end type ScratchFileBufferType + +contains + + subroutine scratch_init(this) + class(ScratchFileBufferType) :: this + integer(I4B) :: istat + open (newunit=this%iun, status='scratch', form='unformatted', & + access='sequential', iostat=istat) + if (istat /= 0) & + call pstop(1, 'failed to open scratch track buffer file') + end subroutine scratch_init + + subroutine scratch_append(this, rec) + class(ScratchFileBufferType) :: this + type(ParticleTrackRecordType), intent(in) :: rec + write (this%iun) rec + this%nrecords = this%nrecords + 1 + end subroutine scratch_append + + subroutine scratch_flush(this, files) + class(ScratchFileBufferType) :: this + type(ParticleTrackFileType), intent(in) :: files(:) + integer(I4B) :: n, i + type(ParticleTrackRecordType) :: rec + + rewind (this%iun) + do n = 1, this%nrecords + read (this%iun) rec + do i = 1, size(files) + if (files(i)%iun > 0 .and. & + (files(i)%iprp == -1 .or. files(i)%iprp == rec%iprp)) & + call save_record(files(i)%iun, rec, files(i)%csv) + end do + end do + rewind (this%iun) + this%nrecords = 0 + end subroutine scratch_flush + + subroutine scratch_discard(this) + class(ScratchFileBufferType) :: this + rewind (this%iun) + this%nrecords = 0 + end subroutine scratch_discard + + subroutine scratch_destroy(this) + class(ScratchFileBufferType) :: this + if (this%iun > 0) then + close (this%iun) + this%iun = 0 + end if + this%nrecords = 0 + end subroutine scratch_destroy + +end module ScratchFileBufferModule diff --git a/src/meson.build b/src/meson.build index a281db6491a..736f4d110ff 100644 --- a/src/meson.build +++ b/src/meson.build @@ -363,6 +363,9 @@ modflow_sources = files( 'Solution' / 'ParticleTracker' / 'Method' / 'TernarySolveUtils.f90', 'Solution' / 'ParticleTracker' / 'Particle' / 'Particle.f90', 'Solution' / 'ParticleTracker' / 'Particle' / 'ParticleReleaseSchedule.f90', + 'Solution' / 'ParticleTracker' / 'Particle' / 'ParticleTrackEventBuffer.f90', + 'Solution' / 'ParticleTracker' / 'Particle' / 'MemoryBuffer.f90', + 'Solution' / 'ParticleTracker' / 'Particle' / 'ScratchFileBuffer.f90', 'Solution' / 'ParticleTracker' / 'Particle' / 'ParticleTracks.f90', 'Solution' / 'ParticleTracker' / 'Particle' / 'Events' / 'ParticleEvents.f90', 'Solution' / 'ParticleTracker' / 'Particle' / 'Events' / 'ParticleEvent.f90',