diff --git a/autotest/test_prt_dry.py b/autotest/test_prt_dry.py index 1b08c19d2ef..6cb0a138946 100644 --- a/autotest/test_prt_dry.py +++ b/autotest/test_prt_dry.py @@ -1,15 +1,13 @@ """ Tests particle tracking in "dry" conditions. -Particles should terminate in inactive cells. - The PRP package provides the `DRY` option to specify how particles should behave in dry conditions when the flow model enables the Newton formulation. This test case is adapted from the example -simulation provided by Javier Gonzalez in +simulation provided by @javgs-bd in https://github.com/MODFLOW-USGS/modflow6/issues/2014. """ diff --git a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn index 19b8210a510..76aa68560c8 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn @@ -173,7 +173,7 @@ type string valid drop stop stay reader urword optional true -longname what to do in active but dry cells +longname what to do in dry-but-active cells description is a string indicating how particles should behave in dry-but-active cells (as can occur with the Newton formulation). The value can be ``DROP'', ``STOP'', or ``STAY''. The default is ``DROP'', which passes particles vertically and instantaneously to the water table. ``STOP'' causes particles to terminate. ``STAY'' causes particles to remain stationary but active. block options diff --git a/doc/mf6io/prt/prt.tex b/doc/mf6io/prt/prt.tex index 372004e7438..385f45f4260 100644 --- a/doc/mf6io/prt/prt.tex +++ b/doc/mf6io/prt/prt.tex @@ -29,26 +29,20 @@ \subsection{Vertical Tracking} Normally, an inactive cell might be dry or explicitly disabled (idomain). With Newton, dry cells remain active. -Of the flow model, PRT knows only what the FMI or exchange tells it. This includes heads, flows and the active grid region, but not whether Newton is on, where boundary packages are, or even which boundary packages are present. - -Tracking and termination decisions must be made, therefore, on the basis of available information like whether a cell is active, whether the cell is dry, and whether the particle is dry (above the water table). - Release-time and tracking-time considerations are described (and implemented) separately. \subsubsection{Release} At release time, PRT decides whether to release each particle or to terminate it unreleased. -If the release cell is active, the particle will be released at the specified coordinates. +If the release cell is active, the particle will be released at the user-specified location. If the release cell is inactive, behavior is determined by the DRAPE option. If the DRAPE option is enabled, the particle will be released from the top-most active cell beneath it, if any. If there is no active cell underneath the particle in any layer, or if DRAPE is not enabled, the particle will terminate unreleased (with status code 8). -Since under the Newton formulation dry cells can remain active, the DRAPE option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured with tracking-time settings. +Since under the Newton formulation dry cells can remain active, the DRAPE option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured using the DRY\_TRACKING\_METHOD option discussed below. \subsubsection{Tracking} -A particle might find itself above the water table for one of two reasons: it was released above the water table, or the water table has receded. - With the Newton formulation, particles can be released into dry-but-active cells. Particle trajectories are solved over the same time discretization used by the flow model. A particle may be immersed in the flow field in one time step, and find that the water table has dropped below it in the next. @@ -57,11 +51,11 @@ \subsubsection{Tracking} A particle in a dry but active cell, or above the water table in a partially saturated cell, need not terminate. We call such a particle dry. The PRP package provides an option DRY\_TRACKING\_METHOD determining how dry particles should behave. Supported values are DROP (default), STOP, and STAY. -If DROP is selected, or if a DRY\_TRACKING\_METHOD is unspecified, a particle in a dry position is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. If the vertical column containing the particle is entirely dry, the particle will terminate upon reaching the bottom of the model grid. +If DROP is selected, or if a DRY\_TRACKING\_METHOD is unspecified, a dry particle is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. If the vertical column containing the particle is entirely dry, the particle will terminate upon reaching the bottom of the model grid. If STOP is selected, dry particles will be terminated. -If STAY is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue or b) the simulation ends. +If STAY is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue, b) the particle terminates due to reaching its STOPTIME or STOPTRAVELTIME, or c) the simulation ends. \subsection{Particle Track Output} diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index 4298fcca863..0730d8f6b02 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/src/Model/ParticleTracking/prt.f90 @@ -957,7 +957,7 @@ subroutine prt_solve(this) ! Get and apply the tracking method call this%method%apply(particle, tmax) - ! Reset previous cell and zone numbers + ! Reset previous cell, exit face, and zone numbers particle%icp = 0 particle%izp = 0 diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index 3700e3e80ce..58648e42f45 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -199,11 +199,11 @@ subroutine load_particle(this, cell, particle) ipos = idiag + inbr ic = dis%con%ja(ipos) - if (ic == particle%icp) then - ! terminate in the previous cell. - ! got here by falling through the - ! bottom of a well.. TODO: check - ! if entry/exit face is top/bot? + ! if returning to immediately previous cell + ! through its bottom, we've entered a cycle + ! as can occur e.g. in the bottom of wells. + ! terminate in the previous cell. + if (ic == particle%icp .and. inface == 7) then particle%advancing = .false. particle%idomain(2) = particle%icp particle%istatus = 2 diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index d33a60e1ba0..b1dbcd54419 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -156,11 +156,11 @@ subroutine load_particle(this, cell, particle) ipos = idiag + inbr ic = dis%con%ja(ipos) - if (ic == particle%icp) then - ! terminate in the previous cell. - ! got here by falling through the - ! bottom of a well.. TODO: check - ! if entry/exit face is top/bot? + ! if returning to immediately previous cell + ! through its bottom, we've entered a cycle + ! as can occur e.g. in the bottom of wells. + ! terminate in the previous cell. + if (ic == particle%icp .and. inface == 7) then particle%advancing = .false. particle%idomain(2) = particle%icp particle%istatus = 2 diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index add4dbac914..347ae4b1a2b 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -311,10 +311,7 @@ end subroutine save_particle !! !! Apply a translation and/or rotation to particle coordinates. !! No rescaling. It's also possible to invert a transformation. - !! - !! Upon inverting, be sure to untransform to eliminate roundoff - !! error accumulated in the particle origin/rotation variables, - !! and unsetting the transformed flag. + !! Be sure to reset the transformation after using it. !< subroutine transform_coords(this, xorigin, yorigin, zorigin, & sinrot, cosrot, invert)