fix(prt): support exchange coupling to ATS flow models#2753
Conversation
d24b5ec to
4278c86
Compare
| [[items]] | ||
| section = "fixes" | ||
| subsection = "model" | ||
| description = "A PRT model could not previously be coupled via an exchange to a GWF model that uses adaptive time stepping (ATS). This is now supported." |
There was a problem hiding this comment.
is the issue that it could not be coupled, or that when coupled to a GWF-ATS model it did not work properly (in that it did not reset the particles)?
There was a problem hiding this comment.
Good point, it could be coupled but the results would be incorrect. Will improve the phrasing.
| do i = 1, this%ntrackfiles | ||
| file = this%files(i) | ||
| if (file%iun > 0 .and. & | ||
| (file%iprp == -1 .or. file%iprp == rec%iprp)) & |
There was a problem hiding this comment.
It looks like the ParticleTracks object is in charge of determining where the buffered records are written. To do this, it needs to be aware of the kinds of things that can have an associated output file -- currently, the model (iprp == -1 ?) and the prp's. just wondering whether there's a (reasonable) way to arrange this so that ParticleTracks doesn't need to be aware of that, and whether that would be desirable. I.e., the model and the prp's (and whatever other objects creating track output there might hypothetically be in the future -- advanced packages?) decide for themselves which records they want in their output. Maybe each even has its own buffer? Could be something to chat about at some point.
There was a problem hiding this comment.
Nice idea.. Currently the model owns a single ParticleTracksType which knows about model- and PRP-scoped output files, which is awkward. Would make more sense for each to manage its own output. It would be really nice if model and package could select events independently. Need to consider how to configure that, since event selection is in OC.
| ! members | ||
| type(PrtFmiType), pointer :: fmi => null() !< flow model interface | ||
| type(ParticleStoreType), pointer :: particles => null() !< particle store | ||
| type(ParticleStoreType), pointer :: particles_old => null() !< old particle state (for ATS) |
There was a problem hiding this comment.
From the point of view of anything within PRT, does it matter why the tracking is being retried (e.g., because of ATS, or because of Picard)? If not, then maybe most comments (like this one) don't need to refer specifically to ATS, except perhaps in one place just to cite examples of why you might need this retry capability. In other words, maybe the comments can be less ATS-specific overall.
There was a problem hiding this comment.
Will remove mention of ATS here. It may not be the only use case for the old state snapshot in the future, in any case.
There was a problem hiding this comment.
it does matter why the advance hook is being called though (detail in comment below)
| ! Reset state variable | ||
| irestore = 0 | ||
| if (iFailedStepRetry > 0) irestore = 1 | ||
| ! The model should advance just once per (step, ATS retry) pair. |
There was a problem hiding this comment.
Mentioning ATS and Picard as specific examples of retries is good here. However, instead of calling it a "(step, ATS retry) pair" up front, you might consider just "(step, retry) pair", and maybe even rename iFailedStepRetry to simply iRetry to keep it generic?
There was a problem hiding this comment.
iFailedStepRetry comes from SimVariablesModule, so it'd need a local variable if we want to clean this up (which I'm not against doing)
Re: ATS vs Picard, I tried to clarify the comment. The PR description didn't really nail this either at first, I'll edit that too. I figure it may help to reserve "retry" to mean "ATS retry" and "rerun" or "repeated invocation" or similar for when advance is called more than once. Retry reruns need to run normally, but not Picard-triggered reruns.
| ! Coincident release times are merged to | ||
| ! a single time by the release scheduler. | ||
|
|
||
| ! Save or restore particle state for ATS |
|
|
||
| ! Note: particle tracking output is handled elsewhere | ||
| ! Flush buffered track events to disk. prt_ot is called once per | ||
| ! converged time step from Mf6FinalizeTimestep, outside the ATS |
There was a problem hiding this comment.
"outside the retry loop"?
There was a problem hiding this comment.
maybe just not say anything here, and let it be explained elsewhere?
| !! only mutable state (positions, status, tracking time, etc) without | ||
| !! immutable state (identity, options) saves space and time, but note | ||
| !! that it assumes this store already has all the immutable state; if | ||
| !! this routine is ever used beyond the context of ATS compatibility, |
There was a problem hiding this comment.
"retry (e.g., ATS) compatibility"?
| !> @brief Buffer a particle event for deferred write. | ||
| !! The buffer can be flushed to disk by flush_buffer, | ||
| !! or discarded by discard_buffer (e.g., on a failed | ||
| !! ATS retry). The buffer starts at 64 events, then |
There was a problem hiding this comment.
this reference to ATS as an example is good i think
| end subroutine flush_buffer | ||
|
|
||
| !> @brief Discard buffered events. | ||
| !! Called from prt_ad at the start of each new ATS retry attempt, |
There was a problem hiding this comment.
"each new tracking retry"?
aprovost-usgs
left a comment
There was a problem hiding this comment.
I have some comments/questions, but don't let them hold up the merge unless you feel there's more discussion needed
|
#2789 should go before this |
Support exchange-coupling PRT to GWF models that use adaptive time stepping (ATS). This would previously produce incorrect results, as PRT advance/solve routines are generally stateful where the other models are idempotent.
Make these routines idempotent by "staging" tracking state during each solve and only "committing" it on successful time steps.
This also corrects PRT's behavior when it's solved in the same solution group with a flow model with
mxiter > 1. Flopy puts models in the same solution group by default, but we didn't hit this before now becausemxiter = 1by default.This required fixes in
mem_reallocate()routines to correctly handle the array shrinking case (e.g. particles are released in a time step that then fails, so need to be "unreleased"). Existing uses must have been limited to expanding.Also required buffering track records in memory during the time step, and reporting them when the model OT hook is called for a successful time step. The buffer size starts at 64, doubles lazily as needed, and never shrinks, just overwritten on retried steps. (The starting size seems like an arbitrary choice; the release point count seemed like a reasonable guess, but there is no guaranteed relationship between that and the number of particles ultimately released, so I just made it a constant).
This is a departure from the pattern up to now, which was to write particle events immediately, a deliberate choice motivated by MP7 suffering OOM errors when consolidating internal scratch files before writing the final output file for large models. But since just one time step is buffered at once the max usage should typically still be lower than MP7.
Checklist of items for pull request
ruffon new and modified python scripts in .doc, autotests, doc, distribution, pymake, and utils subdirectories.fprettify