From 21ae56fa6205c3b99d4d89c5746fba76a32d0a6e Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Mon, 29 Jan 2024 13:07:54 -0500 Subject: [PATCH] user-specified tracking times --- autotest/prt/prt_test_utils.py | 8 +- .../{test_prt_disv.py => test_prt_disv1.py} | 112 ++++++++++-------- autotest/prt/test_prt_track_events.py | 16 ++- ..._prt_voronoi01.py => test_prt_voronoi1.py} | 31 +++-- ..._prt_voronoi02.py => test_prt_voronoi2.py} | 2 +- autotest/test_gwf_ats02.py | 2 +- doc/mf6io/mf6ivar/dfn/prt-oc.dfn | 37 +++++- doc/mf6io/mf6ivar/md/mf6ivar.md | 21 ++-- doc/mf6io/mf6ivar/tex/prt-oc-desc.tex | 18 ++- doc/mf6io/mf6ivar/tex/prt-oc-options.dat | 3 +- doc/mf6io/mf6ivar/tex/prt-prp-desc.tex | 4 +- doc/mf6io/mf6ivar/tex/prt-prp-options.dat | 2 +- src/Model/ModelUtilities/TrackData.f90 | 13 +- src/Model/ParticleTracking/prt1.f90 | 9 +- src/Model/ParticleTracking/prt1oc1.f90 | 35 +++++- src/Solution/ParticleTracker/Method.f90 | 5 +- .../ParticleTracker/MethodCellPollock.f90 | 5 +- .../ParticleTracker/MethodCellPollockQuad.f90 | 5 +- .../ParticleTracker/MethodCellTernary.f90 | 5 +- src/Solution/ParticleTracker/MethodDis.f90 | 10 +- src/Solution/ParticleTracker/MethodDisv.f90 | 20 +++- .../ParticleTracker/MethodSubcellPollock.f90 | 51 +++++--- .../ParticleTracker/MethodSubcellTernary.f90 | 47 +++++++- src/Utilities/BlockParser.f90 | 27 ++++- 24 files changed, 351 insertions(+), 137 deletions(-) rename autotest/prt/{test_prt_disv.py => test_prt_disv1.py} (83%) rename autotest/prt/{test_prt_voronoi01.py => test_prt_voronoi1.py} (93%) rename autotest/prt/{test_prt_voronoi02.py => test_prt_voronoi2.py} (99%) diff --git a/autotest/prt/prt_test_utils.py b/autotest/prt/prt_test_utils.py index 9fa75bc9e4b..2078fc70d6b 100644 --- a/autotest/prt/prt_test_utils.py +++ b/autotest/prt/prt_test_utils.py @@ -269,17 +269,17 @@ def plot_nodes_and_vertices( # plot vertices verts = mg.verts - ax.plot(verts[:, 0], verts[:, 1], "bo") + ax.plot(verts[:, 0], verts[:, 1], "bo", alpha=0.5) for i in range(ncpl): x, y = verts[i, 0], verts[i, 1] - ax.annotate(str(i + 1), verts[i, :], color="b") + ax.annotate(str(i + 1), verts[i, :], color="b", alpha=0.5) # plot nodes xc, yc = mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0) for i in range(ncpl): x, y = xc[i], yc[i] - ax.plot(x, y, "o", color="grey") - ax.annotate(str(i + 1), (x, y), color="grey") + ax.plot(x, y, "o", color="grey", alpha=0.5) + ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) # create legend ax.legend( diff --git a/autotest/prt/test_prt_disv.py b/autotest/prt/test_prt_disv1.py similarity index 83% rename from autotest/prt/test_prt_disv.py rename to autotest/prt/test_prt_disv1.py index d270b3a4b24..8978107d763 100644 --- a/autotest/prt/test_prt_disv.py +++ b/autotest/prt/test_prt_disv1.py @@ -1,5 +1,6 @@ """ -Tests particle tracking on a vertex (DISV) grid. +Tests particle tracking on a vertex (DISV) grid +that reduces to a regular grid. Two cases are provided, one with valid release position and cell correspondences, and another @@ -30,8 +31,8 @@ from framework import TestFramework -simname = "prtfmi06" -cases = [f"{simname}", f"{simname}bprp"] +simname = "prtdisv1" +cases = [f"{simname}", f"{simname}bprp", f"{simname}trts"] # model info nlay = 1 @@ -51,6 +52,7 @@ nouter, ninner = 100, 300 hclose, rclose, relax = 1e-9, 1e-3, 0.97 porosity = 0.1 +tracktimes = list(np.linspace(0, 11, 1000)) # vertex grid properties disvkwargs = get_disv_kwargs( @@ -211,6 +213,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6): pname="oc", track_filerecord=[prt_track_file], trackcsv_filerecord=[prt_track_csv_file], + tracktimerecord=tracktimes if "trts" in name else None, ) # create the flow model interface @@ -329,6 +332,8 @@ def check_output(idx, test): # load mf6 pathline results mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + if "trts" in name: + assert len(mf6_pls) == 9189 # make sure pathline df has "name" (boundname) column and default values assert "name" in mf6_pls @@ -359,52 +364,56 @@ def check_output(idx, test): spdis = bud.get_data(text="DATA-SPDIS")[0] qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) - # setup plot - fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10)) - for a in ax: - a.set_aspect("equal") - - # plot mf6 pathlines in map view - pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) - pmv.plot_grid() - pmv.plot_array(hds[0], alpha=0.2) - pmv.plot_vector(qx, qy, normalize=True, color="white") - # set zoom area - # xmin, xmax = 2050, 4800 - # ymin, ymax = 5200, 7550 - plot_nodes_and_vertices(gwf, mg, None, mg.ncpl, ax[0]) - mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"]) - for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): - pl.plot( - title="MF6 pathlines", - kind="line", - x="x", - y="y", - ax=ax[0], - legend=False, - color=cm.plasma(ipl / len(mf6_plines)), - ) - - # plot mp7 pathlines in map view - pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1]) - pmv.plot_grid() - pmv.plot_array(hds[0], alpha=0.2) - pmv.plot_vector(qx, qy, normalize=True, color="white") - mp7_plines = mp7_pls.groupby(["particleid"]) - for ipl, (pid, pl) in enumerate(mp7_plines): - pl.plot( - title="MP7 pathlines", - kind="line", - x="x", - y="y", - ax=ax[1], - legend=False, - color=cm.plasma(ipl / len(mp7_plines)), - ) - - # view/save plot - # plt.show() - plt.savefig(gwf_ws / f"test_{simname}.png") + if "bprp" not in name: + # setup plot + fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10)) + for a in ax: + a.set_aspect("equal") + + # plot mf6 pathlines in map view + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) + pmv.plot_grid() + pmv.plot_array(hds[0], alpha=0.2) + pmv.plot_vector(qx, qy, normalize=True, color="white") + # set zoom area + # xmin, xmax = 2050, 4800 + # ymin, ymax = 5200, 7550 + # plot labeled nodes and vertices + plot_nodes_and_vertices(gwf, mg, None, mg.ncpl, ax[0]) + mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"]) + for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): + pl.plot( + title="MF6 pathlines", + linestyle="None" if "trst" in name else "--", + marker="o", + markersize=2, + x="x", + y="y", + ax=ax[0], + legend=False, + color=cm.plasma(ipl / len(mf6_plines)), + ) + + # plot mp7 pathlines in map view + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1]) + pmv.plot_grid() + pmv.plot_array(hds[0], alpha=0.2) + pmv.plot_vector(qx, qy, normalize=True, color="white") + mp7_plines = mp7_pls.groupby(["particleid"]) + for ipl, (pid, pl) in enumerate(mp7_plines): + pl.plot( + title="MP7 pathlines", + kind="line", + x="x", + y="y", + ax=ax[1], + legend=False, + color=cm.plasma(ipl / len(mp7_plines)), + ) + + # view/save plot + # plt.show() + plt.savefig(gwf_ws / f"test_{simname}.png") # convert mf6 pathlines to mp7 format mf6_pls = to_mp7_pathlines(mf6_pls) @@ -428,8 +437,9 @@ def check_output(idx, test): del mp7_pls["node"] # compare mf6 / mp7 pathline data - assert mf6_pls.shape == mp7_pls.shape - assert np.allclose(mf6_pls, mp7_pls, atol=1e-3) + if "trts" not in name: + assert mf6_pls.shape == mp7_pls.shape + assert np.allclose(mf6_pls, mp7_pls, atol=1e-3) @pytest.mark.parametrize("idx, name", enumerate(cases)) diff --git a/autotest/prt/test_prt_track_events.py b/autotest/prt/test_prt_track_events.py index d32aeeda32b..3acffcd0af8 100644 --- a/autotest/prt/test_prt_track_events.py +++ b/autotest/prt/test_prt_track_events.py @@ -44,7 +44,7 @@ from framework import TestFramework -simname = "prtfmi02" +simname = "prtevnt" cases = [ f"{simname}all", f"{simname}rel", @@ -53,6 +53,7 @@ f"{simname}term", f"{simname}wksk", f"{simname}mult", + f"{simname}trts", ] releasepts_prt = { "a": [ @@ -78,6 +79,7 @@ for i in range(5) ], } +tracktimes = list(np.linspace(0, 50, 1000)) def get_output_events(case_name) -> List[str]: @@ -96,6 +98,8 @@ def get_output_events(case_name) -> List[str]: if "term" in case_name else ["RELEASE", "TERMINATE"] if "mult" in case_name + else ["TRACKTIME"] + if "trts" in case_name else ["ALL"] # default ) @@ -170,6 +174,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): track_filerecord=[prt_track_file], trackcsv_filerecord=[prt_track_csv_file], trackeventrecord=events, + tracktimerecord=tracktimes if "trts" in name else None, ) # create the flow model interface @@ -314,12 +319,11 @@ def check_output(idx, test): if "TIMESTEP" in events: # todo: timestep test case pass - if "TIMESERIES" in events: - # todo: timeseries test case - pass if "WEAKSINK" in events: # todo: weak sink test case pass + if "TRACKTIME" in events: + expected_len += 5324 # hardcoded assert len(mf6_pls) == expected_len # make sure mf6 pathline data have correct @@ -373,7 +377,9 @@ def all_equal(col, val): for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): pl.plot( title="MF6 pathlines", - kind="line", + marker="o", + markersize=2, + linestyle="None", x="x", y="y", ax=ax[0], diff --git a/autotest/prt/test_prt_voronoi01.py b/autotest/prt/test_prt_voronoi1.py similarity index 93% rename from autotest/prt/test_prt_voronoi01.py rename to autotest/prt/test_prt_voronoi1.py index f31df4e7dc4..76d1fb9edb5 100644 --- a/autotest/prt/test_prt_voronoi01.py +++ b/autotest/prt/test_prt_voronoi1.py @@ -33,8 +33,10 @@ from framework import TestFramework -simname = "prtvor01" +simname = "prtvor1" cases = [f"{simname}l2r", f"{simname}welp", f"{simname}weli"] +times = [True, False, False] +tracktimes = list(np.linspace(0, 40000, 100)) xmin = 0.0 xmax = 2000.0 ymin = 0.0 @@ -165,7 +167,7 @@ def build_gwf_sim(name, ws, targets): return sim -def build_prt_sim(name, gwf_ws, prt_ws, targets): +def build_prt_sim(idx, name, gwf_ws, prt_ws, targets): prt_ws = Path(prt_ws) gwf_name = get_model_name(name, "gwf") prt_name = get_model_name(name, "prt") @@ -209,7 +211,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets): prpdata = [ # index, (layer, cell index), x, y, z (i, (0, vgrid.intersect(p[0], p[1])), p[0], p[1], p[2]) - for i, p in enumerate(rpts[1:2]) # first release point crashes + for i, p in enumerate(rpts[1:]) # first release point crashes ] prp_track_file = f"{prt_name}.prp.trk" prp_track_csv_file = f"{prt_name}.prp.trk.csv" @@ -217,7 +219,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets): prt, pname="prp1", filename=f"{prt_name}_1.prp", - nreleasepts=1, # len(prpdata), + nreleasepts=len(prpdata), packagedata=prpdata, perioddata={0: ["FIRST"]}, track_filerecord=[prp_track_file], @@ -232,6 +234,8 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets): pname="oc", track_filerecord=[prt_track_file], trackcsv_filerecord=[prt_track_csv_file], + trackeventrecord=["TRACKTIME"] if times[idx] else ["ALL"], + tracktimerecord=tracktimes if times[idx] else None, ) gwf_budget_file = gwf_ws / f"{gwf_name}.bud" gwf_head_file = gwf_ws / f"{gwf_name}.hds" @@ -254,7 +258,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets): def build_models(idx, test): 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 + idx, test.name, test.workspace, test.workspace / "prt", test.targets ) return gwf_sim, prt_sim @@ -282,17 +286,15 @@ def check_output(idx, test): ) if "l2r" in name: - assert pls.shape == (212, 16) + # assert pls.shape == (212, 16) assert (pls.z == 0.5).all() # no z change # path should be horizontal from left to right assert isclose(min(pls.x), 20, rel_tol=1e-4) - assert isclose(max(pls.x), 1979, rel_tol=1e-4) - y = 21 - assert isclose(min(pls.y), y, rel_tol=1e-4) - assert isclose(max(pls.y), y, rel_tol=1e-4) - assert set(endpts.icell) == {921} + assert isclose(max(pls.x), 1980.571, rel_tol=1e-4) + assert isclose(min(pls.y), 21, rel_tol=1e-4) + assert isclose(max(pls.y), 981, rel_tol=1e-4) - plot_2d = True + plot_2d = False if plot_2d: # plot in 2d with mpl fig = plt.figure(figsize=(16, 10)) @@ -347,13 +349,16 @@ def check_output(idx, test): pl.plot( title=title, kind="line", + linestyle="--", + marker="o", + markersize=2, x="x", y="y", ax=ax, legend=False, color="black", ) - # plt.show() + plt.show() plt.savefig(prt_ws / f"{name}.png") plot_3d = False diff --git a/autotest/prt/test_prt_voronoi02.py b/autotest/prt/test_prt_voronoi2.py similarity index 99% rename from autotest/prt/test_prt_voronoi02.py rename to autotest/prt/test_prt_voronoi2.py index d863e426f62..89678169884 100644 --- a/autotest/prt/test_prt_voronoi02.py +++ b/autotest/prt/test_prt_voronoi2.py @@ -32,7 +32,7 @@ from framework import TestFramework -simname = "prtvor02" +simname = "prtvor2" cases = [simname] xmin = 0.0 xmax = 2000.0 diff --git a/autotest/test_gwf_ats02.py b/autotest/test_gwf_ats02.py index 7e03230de24..7514eeafed5 100644 --- a/autotest/test_gwf_ats02.py +++ b/autotest/test_gwf_ats02.py @@ -202,7 +202,7 @@ def make_plot(test): label=f"Layer {ilay + 1}", ) plt.legend() - plt.show() + # plt.show() def check_output(idx, test): diff --git a/doc/mf6io/mf6ivar/dfn/prt-oc.dfn b/doc/mf6io/mf6ivar/dfn/prt-oc.dfn index 8a8622353e0..f98323693d2 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-oc.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-oc.dfn @@ -269,7 +269,7 @@ in_record true optional false tagged true shape -longname tracking events +longname description keyword indicating that particle events will follow block options @@ -282,7 +282,40 @@ optional false tagged false repeating true longname tracking events -description particle event(s) to track. RELEASE selects particle release events. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle termination events. WEAKSINK selects occasions in which a particle exits a weak sink, a cell which removes some but not all of the inflow from adjacent cells. +description event(s) to track. ALL selects all events. RELEASE selects particle release events. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle termination events. WEAKSINK selects occasions in which a particle exits a weak sink, a cell which removes some but not all of the inflow from adjacent cells. TRACKTIME selects times provided as double precision values to the TRACKTIME option. + +block options +name tracktimerecord +type record tracktime times +shape +reader urword +tagged true +optional true +longname +description + +block options +name tracktime +type keyword +reader urword +in_record true +optional false +tagged true +shape +longname +description keyword indicating tracking times will follow + +block options +name times +type double precision +shape (unknown) +reader urword +in_record true +optional false +tagged false +repeating true +longname tracking times +description times to track, relative to the beginning of the simulation. # --------------------- prt oc period --------------------- diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 54c3825a900..60586a6bef2 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -450,7 +450,7 @@ | GWF | RCHA | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that recharge flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | | GWF | RCHA | OPTIONS | TAS6 | KEYWORD | keyword to specify that record corresponds to a time-array-series file. | | GWF | RCHA | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | -| GWF | RCHA | OPTIONS | TAS6_FILENAME | STRING | defines a time-array-series file defining a time-array series that can be used to assign time-varying values. See the Time-Variable Input section for instructions on using the time-array series capability. | +| GWF | RCHA | OPTIONS | TAS6_FILENAME | STRING | defines a time-array-series file defining a time-array series that can be used to assign time-varying values. See the ``Time-Variable Input'' section for instructions on using the time-array series capability. | | GWF | RCHA | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. | | GWF | RCHA | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the Recharge package. See the ``Observation utility'' section for instructions for preparing observation input files. Tables \ref{table:gwf-obstypetable} and \ref{table:gwt-obstypetable} lists observation type(s) supported by the Recharge package. | | GWF | RCHA | PERIOD | IPER | INTEGER | integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. | @@ -539,7 +539,7 @@ | GWF | MAW | CONNECTIONDATA | ICON | INTEGER | integer value that defines the GWF connection number for this multi-aquifer well connection entry. ICONN must be greater than zero and less than or equal to NGWFNODES for multi-aquifer well IFNO. | | GWF | MAW | CONNECTIONDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. One or more screened intervals can be connected to the same CELLID if CONDEQN for a well is MEAN. The program will terminate with an error if MAW wells using SPECIFIED, THIEM, SKIN, or CUMULATIVE conductance equations have more than one connection to the same CELLID. | | GWF | MAW | CONNECTIONDATA | SCRN_TOP | DOUBLE PRECISION | value that defines the top elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_TOP can be any value and is set to the top of the cell. If CONDEQN is MEAN, SCRN\_TOP is set to the multi-aquifer well connection cell top if the specified value is greater than the cell top. The program will terminate with an error if the screen top is less than the screen bottom. | -| GWF | MAW | CONNECTIONDATA | SCRN_BOT | DOUBLE PRECISION | value that defines the bottom elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_BOT can be any value is set to the bottom of the cell. If CONDEQN is MEAN, SCRN\_BOT is set to the multi-aquifer well connection cell bottom if the specified value is less than the cell bottom. The program will terminate with an error if the screen bottom is greater than the screen top. | +| GWF | MAW | CONNECTIONDATA | SCRN_BOT | DOUBLE PRECISION | value that defines the bottom elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_BOT can be any value and is set to the bottom of the cell. If CONDEQN is MEAN, SCRN\_BOT is set to the multi-aquifer well connection cell bottom if the specified value is less than the cell bottom. The program will terminate with an error if the screen bottom is greater than the screen top. | | GWF | MAW | CONNECTIONDATA | HK_SKIN | DOUBLE PRECISION | value that defines the skin (filter pack) hydraulic conductivity (if CONDEQN for the multi-aquifer well is SKIN, CUMULATIVE, or MEAN) or conductance (if CONDEQN for the multi-aquifer well is SPECIFIED) for each GWF node connected to the multi-aquifer well (NGWFNODES). If CONDEQN is SPECIFIED, HK\_SKIN must be greater than or equal to zero. HK\_SKIN can be any value if CONDEQN is THIEM. Otherwise, HK\_SKIN must be greater than zero. If CONDEQN is SKIN, the contrast between the cell transmissivity (the product of geometric mean horizontal hydraulic conductivity and the cell thickness) and the well transmissivity (the product of HK\_SKIN and the screen thicknesses) must be greater than one in node CELLID or the program will terminate with an error condition; if an error condition occurs, it is suggested that the HK\_SKIN be reduced to a value less than K11 and K22 in node CELLID or the THEIM or MEAN conductance equations be used for these multi-aquifer wells. | | GWF | MAW | CONNECTIONDATA | RADIUS_SKIN | DOUBLE PRECISION | real value that defines the skin radius (filter pack radius) for the multi-aquifer well. RADIUS\_SKIN can be any value if CONDEQN is SPECIFIED or THIEM. If CONDEQN is SKIN, CUMULATIVE, or MEAN, the program will terminate with an error if RADIUS\_SKIN is less than or equal to the RADIUS for the multi-aquifer well. | | GWF | MAW | PERIOD | IPER | INTEGER | integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. | @@ -1225,7 +1225,7 @@ | PRT | PRP | OPTIONS | STOP_AT_WEAK_SINK | KEYWORD | is a text keyword to indicate that a particle is to terminate when it enters a cell that is a weak sink. By default, particles are allowed to pass though cells that are weak sinks. | | PRT | PRP | OPTIONS | ISTOPZONE | INTEGER | integer value defining the stop zone number. If cells have been assigned IZONE values in the GRIDDATA block, a particle terminates if it enters a cell whose IZONE value matches ISTOPZONE. An ISTOPZONE value of zero indicates that there is no stop zone. The default value is zero. | | PRT | PRP | OPTIONS | DRAPE | KEYWORD | 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. | -| PRT | PRP | OPTIONS | REFERENCETIME | DOUBLE PRECISION | real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. | +| PRT | PRP | OPTIONS | RELEASETIME | DOUBLE PRECISION | 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 | NRELEASEPTS | INTEGER | is the number of particle release points. | | PRT | PRP | PACKAGEDATA | IRPTNO | INTEGER | integer value that defines the PRP release point number associated with the specified PACKAGEDATA data on the line. IRPTNO must be greater than zero and less than or equal to NRELEASEPTS. The program will terminate with an error if information for a PRP release point number is specified more than once. | | PRT | PRP | PACKAGEDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. | @@ -1234,7 +1234,7 @@ | PRT | PRP | PACKAGEDATA | ZRPT | DOUBLE PRECISION | real value that defines the z coordinate of the release point in model coordinates. The (x, y, z) location specified for the release point must lie within the cell that corresponds to the specified cellid. | | PRT | PRP | PACKAGEDATA | BOUNDNAME | STRING | name of the release-point cell. BOUNDNAME is an ASCII character variable that can contain as many as 40 characters. If BOUNDNAME contains spaces in it, then the entire name must be enclosed within single quotes. | | PRT | PRP | PERIOD | IPER | INTEGER | integer value specifying the stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block applies only to that stress period. | -| PRT | PRP | PERIOD | RELEASESETTING | KEYSTRING | 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. | +| PRT | PRP | PERIOD | RELEASESETTING | KEYSTRING | 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. | | PRT | PRP | PERIOD | ALL | KEYWORD | keyword to indicate release of particles at the start of all time steps in the period. | | PRT | PRP | PERIOD | FIRST | KEYWORD | keyword to indicate release of particles at the start of the first time step in the period. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps. | | PRT | PRP | PERIOD | FREQUENCY | INTEGER | release particles at the specified time step frequency. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps. | @@ -1248,15 +1248,18 @@ | PRT | OC | OPTIONS | CONCENTRATION | KEYWORD | keyword to specify that record corresponds to concentration. | | PRT | OC | OPTIONS | CONCENTRATIONFILE | STRING | name of the output file to write conc information. | | PRT | OC | OPTIONS | PRINT_FORMAT | KEYWORD | keyword to specify format for printing to the listing file. | -| PRT | OC | OPTIONS | TRACKEVENT | STRING | particle tracking event(s) to include in track output files. Can be ALL, RELEASE, TRANSIT, TIMESTEP, TERMINATE, or WEAKSINK. RELEASE selects particle releases. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle terminations. WEAKSINK selects particle exits from weak sink cells. Events may coincide with other events. | -| PRT | OC | OPTIONS | TRACK | KEYWORD | keyword to specify that record corresponds to track. | -| PRT | OC | OPTIONS | TRACKFILE | STRING | name of the output file to write tracking information. | -| PRT | OC | OPTIONS | TRACKCSV | KEYWORD | keyword to specify that record corresponds to the track CSV. | -| PRT | OC | OPTIONS | TRACKCSVFILE | STRING | name of the comma-separated value (CSV) file to write tracking information. | | PRT | OC | OPTIONS | COLUMNS | INTEGER | number of columns for writing data. | | PRT | OC | OPTIONS | WIDTH | INTEGER | width for writing each number. | | PRT | OC | OPTIONS | DIGITS | INTEGER | number of digits to use for writing a number. | | PRT | OC | OPTIONS | FORMAT | STRING | write format can be EXPONENTIAL, FIXED, GENERAL, or SCIENTIFIC. | +| PRT | OC | OPTIONS | TRACK | KEYWORD | keyword to specify that record corresponds to pathlines. | +| PRT | OC | OPTIONS | TRACKFILE | STRING | name of the output file to write tracking information. | +| PRT | OC | OPTIONS | TRACKCSV | KEYWORD | keyword to specify that record corresponds to the track CSV. | +| PRT | OC | OPTIONS | TRACKCSVFILE | STRING | name of the comma-separated value (CSV) file to write tracking information. | +| PRT | OC | OPTIONS | TRACKEVENT | KEYWORD | keyword indicating that particle events will follow | +| PRT | OC | OPTIONS | EVENTS | STRING (UNKNOWN) | event(s) to track. ALL selects all events. RELEASE selects particle release events. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle termination events. WEAKSINK selects occasions in which a particle exits a weak sink, a cell which removes some but not all of the inflow from adjacent cells. TRACKTIME selects times provided as double precision values to the TRACKTIME option. | +| PRT | OC | OPTIONS | TRACKTIME | KEYWORD | keyword indicating tracking times will follow | +| PRT | OC | OPTIONS | TIMES | DOUBLE PRECISION (UNKNOWN) | times to track, relative to the beginning of the simulation. | | PRT | OC | PERIOD | IPER | INTEGER | integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. | | PRT | OC | PERIOD | SAVE | KEYWORD | keyword to indicate that information will be saved this stress period. | | PRT | OC | PERIOD | PRINT | KEYWORD | keyword to indicate that information will be printed this stress period. | diff --git a/doc/mf6io/mf6ivar/tex/prt-oc-desc.tex b/doc/mf6io/mf6ivar/tex/prt-oc-desc.tex index 0e1e650bccb..12e65aa351a 100644 --- a/doc/mf6io/mf6ivar/tex/prt-oc-desc.tex +++ b/doc/mf6io/mf6ivar/tex/prt-oc-desc.tex @@ -19,9 +19,15 @@ \item \texttt{PRINT\_FORMAT}---keyword to specify format for printing to the listing file. -\item \texttt{trackevent}---particle tracking event(s) to include in track output files. Can be ALL, RELEASE, TRANSIT, TIMESTEP, TERMINATE, or WEAKSINK. RELEASE selects particle releases. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle terminations. WEAKSINK selects particle exits from weak sink cells. Events may coincide with other events. +\item \texttt{columns}---number of columns for writing data. + +\item \texttt{width}---width for writing each number. -\item \texttt{TRACK}---keyword to specify that record corresponds to track. +\item \texttt{digits}---number of digits to use for writing a number. + +\item \texttt{format}---write format can be EXPONENTIAL, FIXED, GENERAL, or SCIENTIFIC. + +\item \texttt{TRACK}---keyword to specify that record corresponds to pathlines. \item \texttt{trackfile}---name of the output file to write tracking information. @@ -29,13 +35,13 @@ \item \texttt{trackcsvfile}---name of the comma-separated value (CSV) file to write tracking information. -\item \texttt{columns}---number of columns for writing data. +\item \texttt{TRACKEVENT}---keyword indicating that particle events will follow -\item \texttt{width}---width for writing each number. +\item \texttt{events}---event(s) to track. ALL selects all events. RELEASE selects particle release events. TRANSIT selects cell-to-cell transitions. TIMESTEP selects transitions between timesteps. TERMINATE selects particle termination events. WEAKSINK selects occasions in which a particle exits a weak sink, a cell which removes some but not all of the inflow from adjacent cells. TRACKTIME selects times provided as double precision values to the TRACKTIME option. -\item \texttt{digits}---number of digits to use for writing a number. +\item \texttt{TRACKTIME}---keyword indicating tracking times will follow -\item \texttt{format}---write format can be EXPONENTIAL, FIXED, GENERAL, or SCIENTIFIC. +\item \texttt{times}---times to track, relative to the beginning of the simulation. \end{description} \item \textbf{Block: PERIOD} diff --git a/doc/mf6io/mf6ivar/tex/prt-oc-options.dat b/doc/mf6io/mf6ivar/tex/prt-oc-options.dat index 66d229f1201..08eea4e8fd6 100644 --- a/doc/mf6io/mf6ivar/tex/prt-oc-options.dat +++ b/doc/mf6io/mf6ivar/tex/prt-oc-options.dat @@ -3,7 +3,8 @@ BEGIN OPTIONS [BUDGETCSV FILEOUT ] [CONCENTRATION FILEOUT ] [CONCENTRATION PRINT_FORMAT COLUMNS WIDTH DIGITS ] - [TRACKEVENT ] [TRACK FILEOUT ] [TRACKCSV FILEOUT ] + [TRACKEVENT ] + [TRACKTIME ] END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex b/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex index ee354612dec..3d0f610e9de 100644 --- a/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex +++ b/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex @@ -25,7 +25,7 @@ \item \texttt{DRAPE}---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. -\item \texttt{referencetime}---real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. +\item \texttt{releasetime}---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. \end{description} \item \textbf{Block: DIMENSIONS} @@ -55,7 +55,7 @@ \begin{description} \item \texttt{iper}---integer value specifying the stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block applies only to that stress period. -\item \texttt{releasesetting}---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. +\item \texttt{releasesetting}---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. \begin{lstlisting}[style=blockdefinition] ALL diff --git a/doc/mf6io/mf6ivar/tex/prt-prp-options.dat b/doc/mf6io/mf6ivar/tex/prt-prp-options.dat index 9def10d2c45..70e8a1a9fa1 100644 --- a/doc/mf6io/mf6ivar/tex/prt-prp-options.dat +++ b/doc/mf6io/mf6ivar/tex/prt-prp-options.dat @@ -7,5 +7,5 @@ BEGIN OPTIONS [STOP_AT_WEAK_SINK] [ISTOPZONE ] [DRAPE] - [REFERENCETIME ] + [RELEASETIME ] END OPTIONS diff --git a/src/Model/ModelUtilities/TrackData.f90 b/src/Model/ModelUtilities/TrackData.f90 index 57ea2cff2cf..6bb256af7fa 100644 --- a/src/Model/ModelUtilities/TrackData.f90 +++ b/src/Model/ModelUtilities/TrackData.f90 @@ -32,6 +32,7 @@ module TrackModule !! 2: TIMESTEP !! 3: TERMINATE !! 4: WEAKSINK + !! 5: TRACKTIME !! !! An arbitrary number of files can be managed, internal !! arrays are resized as needed. @@ -45,6 +46,7 @@ module TrackModule logical(LGP), public :: itracktimestep !< track timestep ends logical(LGP), public :: itrackterminate !< track termination events logical(LGP), public :: itrackweaksink !< track weak sink exit events + logical(LGP), public :: itracktimes !< track user-selected times contains procedure :: expand procedure, public :: init_track_file @@ -88,6 +90,7 @@ module TrackModule ! 2: current time step ended**** ! 3: particle terminated ! 4: particle exited weak sink + ! 5: user-specified tracking time ! ! Each record has an "istatus" property, which is the tracking status; ! e.g., awaiting release, active, terminated. A particle may terminate @@ -260,7 +263,8 @@ subroutine save(this, particle, kper, kstp, reason, level) (this%itracktransit .and. reason == 1) .or. & (this%itracktimestep .and. reason == 2) .or. & (this%itrackterminate .and. reason == 3) .or. & - (this%itrackweaksink .and. reason == 4))) & + (this%itrackweaksink .and. reason == 4) .or. & + (this%itracktimes .and. reason == 5))) & return ! -- For now, only allow reporting from outside the tracking @@ -294,14 +298,17 @@ subroutine set_track_events(this, & transit, & timestep, & terminate, & - weaksink) + weaksink, & + tracktimes) class(TrackControlType) :: this - logical(LGP), intent(in) :: release, transit, timestep, terminate, weaksink + logical(LGP), intent(in) :: release, transit, timestep, & + terminate, weaksink, tracktimes this%itrackrelease = release this%itracktransit = transit this%itracktimestep = timestep this%itrackterminate = terminate this%itrackweaksink = weaksink + this%itracktimes = tracktimes end subroutine set_track_events end module TrackModule diff --git a/src/Model/ParticleTracking/prt1.f90 b/src/Model/ParticleTracking/prt1.f90 index eb762052f29..e0bc6557bb7 100644 --- a/src/Model/ParticleTracking/prt1.f90 +++ b/src/Model/ParticleTracking/prt1.f90 @@ -268,7 +268,8 @@ subroutine prt_ar(this) izone=this%mip%izone, & flowja=this%flowja, & porosity=this%mip%porosity, & - retfactor=this%mip%retfactor) + retfactor=this%mip%retfactor, & + tracktimes=this%oc%tracktimes) this%method => method_dis type is (GwfDisvType) call method_disv%init( & @@ -277,7 +278,8 @@ subroutine prt_ar(this) izone=this%mip%izone, & flowja=this%flowja, & porosity=this%mip%porosity, & - retfactor=this%mip%retfactor) + retfactor=this%mip%retfactor, & + tracktimes=this%oc%tracktimes) this%method => method_disv end select @@ -291,7 +293,8 @@ subroutine prt_ar(this) this%oc%itrktrs, & this%oc%itrktst, & this%oc%itrkter, & - this%oc%itrkwsk) + this%oc%itrkwsk, & + this%oc%itrktrt) end subroutine prt_ar !> @brief Read and prepare (calls package read and prepare routines) diff --git a/src/Model/ParticleTracking/prt1oc1.f90 b/src/Model/ParticleTracking/prt1oc1.f90 index afe6d2e6c62..d75ce1885f6 100644 --- a/src/Model/ParticleTracking/prt1oc1.f90 +++ b/src/Model/ParticleTracking/prt1oc1.f90 @@ -6,8 +6,9 @@ module PrtOcModule use OutputControlModule, only: OutputControlType use OutputControlDataModule, only: OutputControlDataType, ocd_cr use SimVariablesModule, only: errmsg, warnmsg - use MemoryManagerModule, only: mem_allocate, mem_deallocate + use MemoryManagerModule, only: mem_allocate, mem_deallocate, mem_reallocate use MemoryHelperModule, only: create_mem_path + use ArrayHandlersModule, only: ExpandArray implicit none private @@ -24,6 +25,8 @@ module PrtOcModule logical(LGP), pointer :: itrktst => null() !< track timestep events logical(LGP), pointer :: itrkter => null() !< track termination events logical(LGP), pointer :: itrkwsk => null() !< track weak sink exit events + logical(LGP), pointer :: itrktrt => null() !< track user-specified times + real(DP), dimension(:), pointer, contiguous :: tracktimes !< user-specified times contains procedure :: oc_ar @@ -76,6 +79,8 @@ subroutine prt_oc_allocate_scalars(this, name_model) call mem_allocate(this%itrktst, 'ITRACKTST', this%memoryPath) call mem_allocate(this%itrkter, 'ITRACKTER', this%memoryPath) call mem_allocate(this%itrkwsk, 'ITRACKWSK', this%memoryPath) + call mem_allocate(this%itrktrt, 'ITRACKTSR', this%memoryPath) + call mem_allocate(this%tracktimes, 0, 'TRACKTIME', this%memoryPath) this%name_model = name_model this%inunit = 0 @@ -91,6 +96,7 @@ subroutine prt_oc_allocate_scalars(this, name_model) this%itrktst = .false. this%itrkter = .false. this%itrkwsk = .false. + this%itrktrt = .false. end subroutine prt_oc_allocate_scalars @@ -157,6 +163,8 @@ subroutine prt_oc_da(this) call mem_deallocate(this%itrktst) call mem_deallocate(this%itrkter) call mem_deallocate(this%itrkwsk) + call mem_deallocate(this%itrktrt) + call mem_deallocate(this%tracktimes) end subroutine prt_oc_da subroutine prt_oc_read_options(this) @@ -176,7 +184,8 @@ subroutine prt_oc_read_options(this) character(len=:), allocatable :: line integer(I4B) :: ierr integer(I4B) :: ipos - logical :: isfound, found, endOfBlock, eventFound + real(DP) :: d + logical(LGP) :: isfound, found, endOfBlock, eventFound, success type(OutputControlDataType), pointer :: ocdobjptr character(len=LINELENGTH) :: trkevent ! -- formats @@ -194,12 +203,12 @@ subroutine prt_oc_read_options(this) ! -- parse options block if detected if (isfound) then write (this%iout, '(/,1x,a,/)') 'PROCESSING OC OPTIONS' + eventFound = .false. do call this%parser%GetNextLine(endOfBlock) if (endOfBlock) exit call this%parser%GetStringCaps(keyword) found = .false. - eventFound = .false. select case (keyword) case ('BUDGETCSV') call this%parser%GetStringCaps(keyword2) @@ -264,6 +273,7 @@ subroutine prt_oc_read_options(this) this%itrktst = .true. this%itrkter = .true. this%itrkwsk = .true. + this%itrktrt = .true. case ('RELEASE') this%itrkrls = .true. case ('TRANSIT') @@ -274,20 +284,35 @@ subroutine prt_oc_read_options(this) this%itrkter = .true. case ('WEAKSINK') this%itrkwsk = .true. + case ('TRACKTIME') + this%itrktrt = .true. case default write (errmsg, '(2a)') & 'Looking for ALL, RELEASE, TRANSIT, TIMESTEP, & - &TERMINATE, or WEAKSINK. Found: ', & + &TERMINATE, WEAKSINK, or TRACKTIME. Found: ', & trim(adjustl(trkevent)) call store_error(errmsg, terminate=.TRUE.) end select eventFound = .true. end do trackeventloop found = .true. + case ('TRACKTIME') + tseventloop: do + success = .false. + call this%parser%TryGetDouble(d, success) + if (.not. success) exit tseventloop + call mem_reallocate(this%tracktimes, size(this%tracktimes) + 1, & + 'TRACKTIME', this%memoryPath) + this%tracktimes(size(this%tracktimes)) = d + end do tseventloop + ! -- todo: check that times are strictly increasing? + this%itrktrt = .true. + found = .true. case default found = .false. end select + ! -- check if we're done if (.not. found) then do ipos = 1, size(this%ocdobj) ocdobjptr => this%ocdobj(ipos) @@ -313,8 +338,10 @@ subroutine prt_oc_read_options(this) this%itrktst = .true. this%itrkter = .true. this%itrkwsk = .true. + this%itrktrt = .true. end if + ! -- logging write (this%iout, '(1x,a)') 'END OF OC OPTIONS' end if end subroutine prt_oc_read_options diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index 86e38c5044b..a277da54563 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -36,6 +36,7 @@ module MethodModule real(DP), dimension(:), pointer, contiguous, public :: flowja => null() !< pointer to intercell flows real(DP), dimension(:), pointer, contiguous, public :: porosity => null() !< pointer to aquifer porosity real(DP), dimension(:), pointer, contiguous, public :: retfactor => null() !< pointer to retardation factor + real(DP), dimension(:), pointer, contiguous, public :: tracktimes => null() contains ! Implemented in all subtypes procedure(apply), deferred :: apply @@ -68,7 +69,7 @@ end subroutine destroy contains subroutine init(this, fmi, cell, subcell, trackctl, & - izone, flowja, porosity, retfactor) + izone, flowja, porosity, retfactor, tracktimes) class(MethodType), intent(inout) :: this type(PrtFmiType), intent(in), pointer, optional :: fmi class(CellType), intent(in), pointer, optional :: cell @@ -78,6 +79,7 @@ subroutine init(this, fmi, cell, subcell, trackctl, & real(DP), intent(in), pointer, optional :: flowja(:) real(DP), intent(in), pointer, optional :: porosity(:) real(DP), intent(in), pointer, optional :: retfactor(:) + real(DP), intent(in), pointer, optional :: tracktimes(:) if (present(fmi)) this%fmi => fmi if (present(cell)) this%cell => cell @@ -87,6 +89,7 @@ subroutine init(this, fmi, cell, subcell, trackctl, & if (present(flowja)) this%flowja => flowja if (present(porosity)) this%porosity => porosity if (present(retfactor)) this%retfactor => retfactor + if (present(tracktimes)) this%tracktimes => tracktimes end subroutine init !> @brief Track particle through subdomains diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index 7e681a88abe..907d36f88a7 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -61,7 +61,10 @@ subroutine load_mcp(this, particle, next_level, submethod) type is (SubcellRectType) call this%load_subcell(particle, subcell) end select - call method_subcell_plck%init(subcell=this%subcell, trackctl=this%trackctl) + call method_subcell_plck%init( & + subcell=this%subcell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_subcell_plck particle%idomain(next_level) = 1 end subroutine load_mcp diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index 432f7efe799..3fe41610c96 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -61,7 +61,10 @@ subroutine load_mcpq(this, particle, next_level, submethod) type is (SubcellRectType) call this%load_subcell(particle, subcell) end select - call method_subcell_plck%init(subcell=this%subcell, trackctl=this%trackctl) + call method_subcell_plck%init( & + subcell=this%subcell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_subcell_plck end subroutine load_mcpq diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index f84e4801427..e2a3cfcaf3e 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -68,7 +68,10 @@ subroutine load_mct(this, particle, next_level, submethod) type is (SubcellTriType) call this%load_subcell(particle, subcell) end select - call method_subcell_tern%init(subcell=this%subcell, trackctl=this%trackctl) + call method_subcell_tern%init( & + subcell=this%subcell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_subcell_tern end subroutine load_mct diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index fde9e82057c..9600bb9c38f 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -79,7 +79,10 @@ subroutine load_dis(this, particle, next_level, submethod) ! -- If cell is active but dry, select and initialize ! -- pass-to-bottom method and set cell method pointer if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? - call method_cell_ptb%init(cell=this%cell, trackctl=this%trackctl) + call method_cell_ptb%init( & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_ptb else ! -- load rectangular cell (todo: refactor into separate routine) @@ -114,7 +117,10 @@ subroutine load_dis(this, particle, next_level, submethod) cell%vz2 = -cell%defn%faceflow(7) * term ! -- Select and initialize Pollock's method and set method pointer - call method_cell_plck%init(cell=this%cell, trackctl=this%trackctl) + call method_cell_plck%init( & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_plck end if end select diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index a9444a109aa..59f78cf270a 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -85,22 +85,34 @@ subroutine load_disv(this, particle, next_level, submethod) if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? ! -- Cell is active but dry, so select and initialize pass-to-bottom ! -- cell method and set cell method pointer - call method_cell_ptb%init(cell=this%cell, trackctl=this%trackctl) + call method_cell_ptb%init( & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_ptb else ! -- Select and initialize cell method and set cell method pointer if (cell%defn%can_be_rect) then call cell_poly_to_rect(cell, rect) base => rect - call method_cell_plck%init(cell=base, trackctl=this%trackctl) + call method_cell_plck%init( & + cell=base, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_plck else if (cell%defn%can_be_quad) then call cell_poly_to_quad(cell, quad) base => quad - call method_cell_quad%init(cell=base, trackctl=this%trackctl) + call method_cell_quad%init( & + cell=base, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_quad else - call method_cell_tern%init(cell=this%cell, trackctl=this%trackctl) + call method_cell_tern%init( & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) submethod => method_cell_tern end if end if diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index e6ce4bc061a..55da23c71e0 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -75,10 +75,10 @@ end subroutine apply_msp !> @brief Track a particle across a rectangular subcell using Pollock's method !! - !! This subroutine consists partly or entirely of code written by - !! David W. Pollock of the USGS for MODPATH 7. The authors of the present - !! code are responsible for its appropriate application in this context - !! and for any modifications or errors. + !! This subroutine consists partly of code written by + !! David W. Pollock of the USGS for MODPATH 7. PRT's + !! authors take responsibility for its application in + !! this context and for any modifications or errors. !< subroutine track_subcell(this, subcell, particle, tmax) ! modules @@ -91,9 +91,9 @@ subroutine track_subcell(this, subcell, particle, tmax) real(DP), intent(in) :: tmax ! local doubleprecision :: vx, dvxdx, vy, dvydy, vz, dvzdz - doubleprecision :: dtexitx, dtexity, dtexitz, dtexit, texit, dt, t + doubleprecision :: dtexitx, dtexity, dtexitz, dtexit, texit, dt, t, t0 doubleprecision :: x, y, z - integer :: statusVX, statusVY, statusVZ + integer :: statusVX, statusVY, statusVZ, i doubleprecision :: initialX, initialY, initialZ integer :: exitFace integer :: reason @@ -113,14 +113,8 @@ subroutine track_subcell(this, subcell, particle, tmax) statusVZ = calculate_dt(subcell%vz1, subcell%vz2, subcell%dz, & initialZ, vz, dvzdz, dtexitz) - ! -- Check for no exit face + ! -- Subcells should never be strong sinks, contact the developer situation if ((statusVX .eq. 3) .and. (statusVY .eq. 3) .and. (statusVZ .eq. 3)) then - ! exitFace = 0 - ! particle%iboundary(3) = exitFace - ! particle%istatus = 5 - ! return - - ! contact the developer situation (for now? always?) print *, "Subcell with no exit face:", & "particle", get_particle_id(particle), & "cell", particle%idomain(2) @@ -161,13 +155,42 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- Compute exit time texit = particle%ttrack + dtexit + t0 = particle%ttrack + + ! -- Solve for user-specified tracking times + if (associated(this%tracktimes)) then + ! todo: start from period/timestep and/or current simulation time + do i = 1, size(this%tracktimes) + t = this%tracktimes(i) + if (t == particle%ttrack) then + call this%trackctl%save(particle, kper=kper, & + kstp=kstp, reason=5) + cycle + else if (t > particle%ttrack .and. t < texit .and. t < tmax) then + dt = t - t0 + x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, & + dt, initialX, subcell%dx, statusVX == 1) + y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, & + dt, initialY, subcell%dy, statusVY == 1) + z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, & + dt, initialZ, subcell%dz, statusVZ == 1) + particle%x = x * subcell%dx + particle%y = y * subcell%dy + particle%z = z * subcell%dz + particle%ttrack = t + particle%istatus = 1 + call this%trackctl%save(particle, kper=kper, & + kstp=kstp, reason=5) + end if + end do + end if if (texit .gt. tmax) then ! -- The computed exit time is greater than the maximum time, so set ! -- final time for particle trajectory equal to maximum time and ! -- calculate particle location at that final time. t = tmax - dt = t - particle%ttrack + dt = t - t0 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, & dt, initialX, subcell%dx, statusVX == 1) y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, & diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 57a684f83e5..3257ea25c1c 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -79,12 +79,12 @@ subroutine track_subcell(this, subcell, particle, tmax) double precision :: rot(2, 2), res(2), loc(2) double precision :: alp, bet, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti double precision :: vzbot, vztop, vzi, vziodz, az, dtexitz - double precision :: dt, t, dtexitxy, texit, x, y, z + double precision :: dt, t, t0, dtexitxy, texit, x, y, z integer :: izstatus, itopbotexit integer :: ntmax, nsave, isolv, itrifaceenter, itrifaceexit double precision :: tol, step, dtexit, alpexit, betexit integer :: ntdebug ! kluge - integer :: reason + integer :: reason, i lbary = .true. ! kluge ntmax = 10000 @@ -191,12 +191,53 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- Compute exit time texit = particle%ttrack + dtexit + t0 = particle%ttrack + + ! -- Solve for user-specified tracking times + if (associated(this%tracktimes)) then + ! todo: start from period/timestep and/or current simulation time + do i = 1, size(this%tracktimes) + t = this%tracktimes(i) + if (t == particle%ttrack) then + call this%trackctl%save(particle, kper=kper, & + kstp=kstp, reason=5) + cycle + else if (t > particle%ttrack .and. t < texit .and. t < tmax) then + dt = t - t0 + call step_analytical(dt, alp, bet) + loc = (/alp, bet/) + if (lbary) loc = skew(loc, (/sxx, sxy, syy/), invert=.true.) + rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot)) + res = matmul(rot, loc) ! rotate vector + x = res(1) + x0 + y = res(2) + y0 + if (izstatus .eq. 2) then ! kluge note: make this into a function + ! -- vz uniformly zero + z = zi + else if (izstatus .eq. 1) then + ! -- vz uniform, nonzero + z = zi + vzi * dt + else + ! -- vz nonuniform + z = zbot + (vzi * dexp(az * dt) - vzbot) / az + end if + ! end if + particle%x = x + particle%y = y + particle%z = z + particle%ttrack = t + particle%istatus = 1 + call this%trackctl%save(particle, kper=kper, & + kstp=kstp, reason=5) + end if + end do + end if if (texit .gt. tmax) then ! -- The computed exit time is greater than the maximum time, so set ! -- final time for particle trajectory equal to maximum time. t = tmax - dt = t - particle%ttrack + dt = t - t0 exitFace = 0 particle%istatus = 1 particle%advancing = .false. diff --git a/src/Utilities/BlockParser.f90 b/src/Utilities/BlockParser.f90 index 887955ff972..d808748f036 100644 --- a/src/Utilities/BlockParser.f90 +++ b/src/Utilities/BlockParser.f90 @@ -6,7 +6,7 @@ !< module BlockParserModule - use KindModule, only: DP, I4B + use KindModule, only: DP, I4B, LGP use DevFeatureModule, only: dev_feature use ConstantsModule, only: LENBIGLINE, LENHUGELINE, LINELENGTH, MAXCHARLEN use InputOutputModule, only: urword, upcase, openfile, & @@ -39,6 +39,7 @@ module BlockParserModule procedure, public :: GetCellid procedure, public :: GetCurrentLine procedure, public :: GetDouble + procedure, public :: TryGetDouble procedure, public :: GetInteger procedure, public :: GetLinesRead procedure, public :: GetNextLine @@ -313,11 +314,29 @@ function GetDouble(this) result(r) if (istart == istop .and. istop == len(this%line)) then call this%ReadScalarError('DOUBLE PRECISION') end if - ! - ! -- return - return + end function GetDouble + subroutine TryGetDouble(this, r, success) + ! -- dummy variables + class(BlockParserType), intent(inout) :: this !< BlockParserType object + real(DP), intent(inout) :: r !< double precision real variable + logical(LGP), intent(inout) :: success !< whether parsing was successful + ! -- local variables + integer(I4B) :: istart + integer(I4B) :: istop + integer(I4B) :: ival + + call urword(this%line, this%lloc, istart, istop, 3, ival, r, & + this%iout, this%iuext) + + success = .true. + if (istart == istop .and. istop == len(this%line)) then + success = .false. + end if + + end subroutine TryGetDouble + !> @ brief Issue a read error !! !! Method to issue an unable to read error.