Skip to content

Commit

Permalink
refactor prp:
Browse files Browse the repository at this point in the history
* simplify implementation
* only support scalar time offset (fraction)
* rename referencetime to releasetime
* check particle release location elevation
* modify array expansion to respect lower bound
  • Loading branch information
wpbonelli committed Dec 18, 2023
1 parent 9993033 commit 8f9b849
Show file tree
Hide file tree
Showing 10 changed files with 414 additions and 632 deletions.
30 changes: 18 additions & 12 deletions autotest/prt/test_prt_fmi05.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
"""
Test cases exercising release timing, 1st via
package-level REFERENCETIME option, then with
package-level RELEASETIME option, & then with
period-block config STEPS 1 and FRACTION 0.5.
The model is setup to release halfway through
the first and only time step of the first and
only stress period, with duration 1 time unit,
so the same value of 0.5 can be used for both
REFERENCETIME and FRACTION.
RELEASETIME and FRACTION.
Period-block FRACTION should work with FIRST
and ALL, but flopy hangs with either option.
Expand All @@ -18,7 +18,7 @@
Particles are released from the top left cell.
Results are compared against a MODPATH 7 model.
Referencetime 0.5 could be configured, but mp7
Telease time 0.5 could be configured, but mp7
reports relative times, so there is no reason
& mp7 results are converted before comparison.
"""
Expand All @@ -37,13 +37,19 @@
from flopy.plot.plotutil import to_mp7_pathlines
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from prt_test_utils import (all_equal, check_budget_data, check_track_data,
get_gwf_sim, get_model_name, get_partdata)
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_gwf_sim,
get_model_name,
get_partdata,
)

simname = "prtfmi05"
ex = [
cases = [
# options block options
f"{simname}reft", # REFERENCETIME 0.5
f"{simname}relt", # RELEASETIME 0.5
# period block options
# f"{simname}all", # ALL FRACTION 0.5 # todo debug flopy hanging
# f"{simname}frst", # FIRST FRACTION 0.5 # todo debug flopy hanging
Expand All @@ -52,7 +58,7 @@


def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]:
if "reft" in name:
if "relt" in name:
return None
opt = [
"FIRST"
Expand Down Expand Up @@ -115,8 +121,8 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
prp_track_file = f"{prtname}.prp.trk"
prp_track_csv_file = f"{prtname}.prp.trk.csv"
pdat = get_perioddata(prtname, fraction=fraction)
# fraction 0.5 equiv. to reftime 0.5 since 1 period 1 step with length 1
reft = fraction if "reft" in prtname else None
# fraction 0.5 equiv. to release time 0.5 since 1 period 1 step with length 1
trelease = fraction if "relt" in prtname else None
flopy.mf6.ModflowPrtprp(
prt,
pname="prp1",
Expand All @@ -126,7 +132,7 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
perioddata=pdat,
track_filerecord=[prp_track_file],
trackcsv_filerecord=[prp_track_csv_file],
referencetime=reft,
releasetime=trelease,
)

# create output control package
Expand Down Expand Up @@ -195,7 +201,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
return mp


@pytest.mark.parametrize("name", ex)
@pytest.mark.parametrize("name", cases)
@pytest.mark.parametrize("fraction", [0.5])
def test_mf6model(name, function_tmpdir, targets, fraction):
# workspace
Expand Down
13 changes: 10 additions & 3 deletions autotest/prt/test_prt_fmi06.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,14 @@
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.utils.gridutil import get_disv_kwargs
from prt_test_utils import (all_equal, check_budget_data, check_track_data,
get_partdata, has_default_boundnames,
plot_nodes_and_vertices)
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_partdata,
has_default_boundnames,
plot_nodes_and_vertices,
)

simname = "prtfmi06"
ex = [f"{simname}", f"{simname}bprp"]
Expand Down Expand Up @@ -442,5 +447,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
del mp7_pls["node"]

# compare mf6 / mp7 pathline data
# import pdb
# pdb.set_trace()
assert mf6_pls.shape == mp7_pls.shape
assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
8 changes: 4 additions & 4 deletions doc/mf6io/mf6ivar/dfn/prt-prp.dfn
Original file line number Diff line number Diff line change
Expand Up @@ -128,12 +128,12 @@ longname drape
description is a text keyword to indicate that if a particle's release point is in a cell that happens to be dry at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is dry at release time, instead the particle is terminated immediately with ireason 3 and istatus 8.

block options
name referencetime
name releasetime
type double precision
reader urword
optional true
longname period referencetime
description real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1.
longname global release time
description real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. Overrides release settings specified in period data.

# --------------------- prt prp dimensions ---------------------

Expand Down Expand Up @@ -248,7 +248,7 @@ tagged false
in_record true
reader urword
longname
description specifies when to release particles within the stress period. Overrides package-level REFERENCETIME option, which applies to all stress periods. By default, RELEASESETTING configures particles for release at the beginning of the specified time steps. For time-offset releases, provide a FRACTION with another RELEASESETTING option.
description specifies when to release particles within the stress period. Overrides package-level RELEASETIME option, which applies to all stress periods. By default, RELEASESETTING configures particles for release at the beginning of the specified time steps. For time-offset releases, provide a FRACTION value.

block period
name all
Expand Down
2 changes: 1 addition & 1 deletion src/Model/ExplicitModel.f90
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ module ExplicitModelModule
!<
type, extends(BaseModelType) :: ExplicitModelType
character(len=LINELENGTH), pointer :: filename => null() !< input file name
integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< ibound
integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< ibound
type(ListType), pointer :: bndlist => null() !< array of boundary packages
class(DisBaseType), pointer :: dis => null() !< discretization object
contains
Expand Down
8 changes: 4 additions & 4 deletions src/Model/ParticleTracking/prt1.f90
Original file line number Diff line number Diff line change
Expand Up @@ -432,7 +432,7 @@ subroutine prt_cq_sto(this)
packobj => GetBndFromList(this%bndlist, ip)
select type (packobj)
type is (PrtPrpType)
do np = 1, packobj%npart
do np = 1, packobj%nparticles
istatus = packobj%particles%istatus(np)
! refine these conditions as necessary
! (status 8 is permanently unreleased)
Expand Down Expand Up @@ -952,9 +952,9 @@ subroutine prt_solve(this)
end if

! -- Loop over particles in package
do np = 1, packobj%npart
do np = 1, packobj%nparticles
! -- Load particle from storage
call particle%load_from_list(packobj%particles, this%id, iprp, np)
call particle%load_from_store(packobj%particles, this%id, iprp, np)

! -- If particle is permanently unreleased, record its initial/terminal state
if (particle%istatus == 8) &
Expand All @@ -980,7 +980,7 @@ subroutine prt_solve(this)
call this%method%apply(particle, tmax)

! -- Update particle storage
call packobj%particles%update_from_particle(particle, np)
call packobj%particles%load_from_particle(particle, np)
end do
end select
end do
Expand Down
Loading

0 comments on commit 8f9b849

Please sign in to comment.