From 16ad3180a9e77c59e9c574f96c21343d038fa891 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Tue, 25 Jun 2024 19:40:47 -0400 Subject: [PATCH] feat(prp): support local, offset and transformed coords --- autotest/test_prt_coord_transform.py | 235 ++++++++++++++++++ autotest/test_prt_fmi.py | 1 - doc/ReleaseNotes/develop.tex | 2 +- doc/mf6io/mf6ivar/dfn/prt-prp.dfn | 26 +- src/Model/ModelUtilities/TrackControl.f90 | 5 +- src/Model/ModelUtilities/TrackFile.f90 | 23 +- src/Model/ParticleTracking/prt-prp.f90 | 83 ++++++- src/Solution/ParticleTracker/Method.f90 | 5 +- .../ParticleTracker/MethodCellPollock.f90 | 1 + .../ParticleTracker/MethodCellPollockQuad.f90 | 2 +- .../ParticleTracker/MethodCellTernary.f90 | 1 - src/Solution/ParticleTracker/Particle.f90 | 60 +++-- 12 files changed, 401 insertions(+), 43 deletions(-) create mode 100644 autotest/test_prt_coord_transform.py diff --git a/autotest/test_prt_coord_transform.py b/autotest/test_prt_coord_transform.py new file mode 100644 index 00000000000..8acf009b325 --- /dev/null +++ b/autotest/test_prt_coord_transform.py @@ -0,0 +1,235 @@ +""" +Tests ability to transform release point coordinates from +global, local, or local offset coordinate representations +to model coordinates, as well as to transform output data +from model to global cooordinates. + +The grid is a 10x10 square with a single layer, +the same flow system shown on the FloPy readme. + +Test cases are defined for each coordinate system option. +These are: + + - local_xy + - local_xy_offset + - dev_global_xy + +""" + +from pathlib import Path + +import flopy +import matplotlib.cm as cm +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +from flopy.utils import PathlineFile +from flopy.utils.binaryfile import HeadFile +from flopy.utils.gridintersect import GridIntersect +from shapely.geometry import Point +from framework import TestFramework +from prt_test_utils import ( + FlopyReadmeCase, + all_equal, + check_budget_data, + check_track_data, + get_model_name, + get_partdata, + has_default_boundnames, + DEFAULT_EXIT_SOLVE_TOL, +) + +simname = "prtcrd" +cases = [f"{simname}mdl", f"{simname}lcl", f"{simname}ofs"] # todo global +nodes = list(range(100)) + + +def get_start_points(): + return np.add( + np.transpose( + np.array(np.meshgrid(range(10), range(10))).reshape(2, -1) + ), + 0.5, + ) + + +def get_prp(prt): + gi = GridIntersect(prt.modelgrid, method="vertex", rtree=True) + if "mdl" in prt.name: + startpts = get_start_points() + releasepts = [ + (i, (0, *gi.intersect(Point(rpt))[0].cellids), *rpt, 0.5) + for i, rpt in enumerate(startpts) + ] + elif "lcl" in prt.name: + releasepts = [ + (nn, tuple([*prt.modelgrid.get_lrc([nn])[0]]), 0.5, 0.5, 0.5) + for nn in nodes + ] + elif "ofs" in prt.name: + releasepts = [ + (nn, tuple([*prt.modelgrid.get_lrc([nn])[0]]), 0.0, 0.0, 0.0) + for nn in nodes + ] + elif "gbl" in prt.name: + # todo global + pass + + prp_track_file = f"{prt.name}_prt.prp.trk" + prp_track_csv_file = f"{prt.name}_prt.prp.trk.csv" + return flopy.mf6.ModflowPrtprp( + prt, + pname="prp1", + filename=f"{prt.name}_1.prp", + nreleasepts=len(releasepts), + packagedata=releasepts, + perioddata={0: ["FIRST"]}, + track_filerecord=[prp_track_file], + trackcsv_filerecord=[prp_track_csv_file], + local_xy="lcl" in prt.name, + local_xy_offset="ofs" in prt.name, + exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, + ) + + +def build_prt_sim(name, gwf_ws, prt_ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=prt_ws, + ) + + # create tdis package + flopy.mf6.modflow.mftdis.ModflowTdis( + sim, + pname="tdis", + time_units="DAYS", + nper=FlopyReadmeCase.nper, + perioddata=[ + ( + FlopyReadmeCase.perlen, + FlopyReadmeCase.nstp, + FlopyReadmeCase.tsmult, + ) + ], + ) + + # create prt model + prt_name = get_model_name(name, "prt") + prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name, save_flows=True) + + # create prt discretization + flopy.mf6.modflow.mfgwfdis.ModflowGwfdis( + prt, + pname="dis", + nlay=FlopyReadmeCase.nlay, + nrow=FlopyReadmeCase.nrow, + ncol=FlopyReadmeCase.ncol, + ) + + # create mip package + flopy.mf6.ModflowPrtmip( + prt, pname="mip", porosity=FlopyReadmeCase.porosity + ) + + # create prp package + prp = get_prp(prt) + + # create output control package + prt_budget_file = f"{prt_name}.bud" + prt_track_file = f"{prt_name}.trk" + prt_track_csv_file = f"{prt_name}.trk.csv" + flopy.mf6.ModflowPrtoc( + prt, + pname="oc", + budget_filerecord=[prt_budget_file], + track_filerecord=[prt_track_file], + trackcsv_filerecord=[prt_track_csv_file], + saverecord=[("BUDGET", "ALL")], + ) + + # create the flow model interface + gwf_name = get_model_name(name, "gwf") + 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), + ], + ) + + # add explicit model solution + ems = flopy.mf6.ModflowEms( + sim, + pname="ems", + filename=f"{prt_name}.ems", + ) + sim.register_solution_package(ems, [prt.name]) + + return sim + + +def build_models(idx, test): + gwf_ws = test.workspace / "gwf" + prt_ws = test.workspace / "prt" + gwf_sim = FlopyReadmeCase.get_gwf_sim( + test.name, gwf_ws, test.targets["mf6"] + ) + prt_sim = build_prt_sim(test.name, gwf_ws, prt_ws, test.targets["mf6"]) + return gwf_sim, prt_sim + + +def check_output(idx, test): + name = test.name + gwf_ws = test.workspace + prt_ws = test.workspace / "prt" + gwf_sim, prt_sim = test.sims[:2] + gwf = gwf_sim.get_model() + prt = prt_sim.get_model() + hds = gwf.output.head().get_data() + mg = gwf.modelgrid + pathlines = pd.read_csv(prt_ws / f"{prt.name}.trk.csv") + startpts = pathlines[pathlines.ireason == 0] + endpts = pathlines[pathlines.ireason == 3] + + # check start points + assert np.allclose( + np.sort(startpts[["x", "y"]].to_numpy(), axis=0), + np.sort(np.array(get_start_points()), axis=0), + ) + + # setup plot + plot_data = False + if plot_data: + fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10)) + for a in ax: + a.set_aspect("equal") + + # plot pathline in map view + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) + pmv.plot_grid() + pmv.plot_array(hds[0], alpha=0.1, cmap="jet") + pmv.plot_pathline(pathlines) + + # view/save plot + plt.show() + plt.savefig(gwf_ws / f"test_{simname}.png") + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + targets=targets, + compare=None, + ) + test.run() diff --git a/autotest/test_prt_fmi.py b/autotest/test_prt_fmi.py index 37ae8d32a71..d2cfdc20f6e 100644 --- a/autotest/test_prt_fmi.py +++ b/autotest/test_prt_fmi.py @@ -1,7 +1,6 @@ """ Tests ability to run a GWF model then a PRT model in separate simulations via flow model interface, -as well as The grid is a 10x10 square with a single layer, the same flow system shown on the FloPy readme. diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 42bd22b4fdc..08e3bb20c1e 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -46,7 +46,7 @@ \item The Stress Package Concentration (SPC) utility available with the SSM Package is now the referred to as the Stress Package Component utility in the MF6IO guide. Additionally, some relatively minor refactoring of the code facilitates use of the SPC utility with the GWE model type so that TEMPERATURE arrays may be read by the utility. The SPC acronym was maintained to preserve backward compatibility. \item The GWF-GWF Exchange has been fixed to support the configuration of the Mover package (MVR) also in cases where the number of exchanges equals zero (NEXG = 0). In earlier versions this has caused MODFLOW to terminate with an error. \item A PRT bug was fixed in which array-based input read from the RCH (recharge) or EVT (evapotranspiration) packages could fail to be processed correctly by the PRT FMI (flow model interface) package, causing a crash. - + \item Previously the PRT model's PRP package required release points to be specified in the model coordinate system (except for the z coordinate via the LOCAL\_Z option). PRT also previously reported pathlines in the model coordinate system. The PRT model now supports coordinate transformations for release point input and pathline output via keyword options for the PRP package, while preserving the existing default model coordinate system. Release point x and y coordinates may now be specified for structured grids with reference to the local rectilinear cell (with coordinates scaled to the unit interval) via the LOCAL\_XY keyword option. The LOCAL\_XY\_OFFSET option, supported for structured or unstructured grids, can be used to locate release point x and y coordinates at an offset from the cell center (with no distance rescaling). \end{itemize} %\underline{INTERNAL FLOW PACKAGES} diff --git a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn index bb5edd48604..ad87ae98163 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn @@ -40,7 +40,31 @@ type keyword reader urword optional true longname whether to use local z coordinates -description indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. +description indicates that ``zrpt'' defines the local z coordinate of the release point within the cell scaled to the unit interval, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. + +block options +name local_xy +type keyword +reader urword +optional true +longname whether to use local x and y coordinates +description indicates that ``xrpt'' and ``yrpt'' define the local x and y coordinates of the release point within the cell scaled to the unit interval, with value of 0 at the west and south faces and 1 at the east and north faces of the cell. This option may only be used with structured models. + +block options +name local_xy_offset +type keyword +reader urword +optional true +longname whether to use local x and y offsets +description indicates that ``xrpt'' and ``rpt'' define the local x and y offsets of the release point within the cell, where the local coordinate system's origin is the cell center and distances are not rescaled. This option may be used with structured or unstructured models. + +block options +name dev_global_xy +type keyword +reader urword +optional true +longname whether to transform x and y coordinates +description indicates that ``xrpt'' and ``yrpt'' define global coordinates and should be transformed to the model coordinate system if georeferencing information (origin and rotation) is available on the model grid. By default, release points are assumed to be in model coordinates. If this option is enabled, particle pathlines will also be reported in global coordinates in output files. block options name extend_tracking diff --git a/src/Model/ModelUtilities/TrackControl.f90 b/src/Model/ModelUtilities/TrackControl.f90 index d16474b471f..7e2e434e3ac 100644 --- a/src/Model/ModelUtilities/TrackControl.f90 +++ b/src/Model/ModelUtilities/TrackControl.f90 @@ -110,10 +110,11 @@ end subroutine expand !! any PRP-level files with PRP index matching the particle's !! PRP index. !< - subroutine save(this, particle, kper, kstp, reason, level) + subroutine save(this, particle, dis, kper, kstp, reason, level) ! dummy class(TrackControlType), intent(inout) :: this type(ParticleType), pointer, intent(in) :: particle + class(DisBaseType), pointer, intent(in) :: dis integer(I4B), intent(in) :: kper integer(I4B), intent(in) :: kstp integer(I4B), intent(in) :: reason @@ -146,7 +147,7 @@ subroutine save(this, particle, kper, kstp, reason, level) if (file%iun > 0 .and. & (file%iprp == -1 .or. & file%iprp == particle%iprp)) & - call save_record(file%iun, particle, & + call save_record(file%iun, particle, dis, & kper, kstp, reason, csv=file%csv) end do end subroutine save diff --git a/src/Model/ModelUtilities/TrackFile.f90 b/src/Model/ModelUtilities/TrackFile.f90 index 4ba94b71e9c..1eaaa8dafa1 100644 --- a/src/Model/ModelUtilities/TrackFile.f90 +++ b/src/Model/ModelUtilities/TrackFile.f90 @@ -2,6 +2,7 @@ module TrackFileModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: DZERO, DPIO180 use ParticleModule, only: ParticleType + use BaseDisModule, only: DisBaseType use GeomUtilModule, only: transform implicit none @@ -85,21 +86,39 @@ module TrackFileModule contains !> @brief Save a particle track record to a binary or CSV file. - subroutine save_record(iun, particle, kper, kstp, reason, csv) + subroutine save_record(iun, particle, dis, kper, kstp, reason, csv) ! dummy integer(I4B), intent(in) :: iun type(ParticleType), pointer, intent(in) :: particle + class(DisBaseType), pointer, intent(in) :: dis integer(I4B), intent(in) :: kper integer(I4B), intent(in) :: kstp integer(I4B), intent(in) :: reason logical(LGP), intent(in) :: csv ! local - real(DP) :: x, y, z + real(DP) :: x, y, z, xout, yout, zout integer(I4B) :: status ! Convert from cell-local to model coordinates if needed call particle%get_model_coords(x, y, z) + ! Apply model grid offset/rotation if needed + if (particle%icoords == 1 .and. & + (dis%angrot /= DZERO .or. & + dis%xorigin /= DZERO .or. & + dis%yorigin /= DZERO)) then + call transform(x, y, z, & + xout, yout, zout, & + dis%xorigin, & + dis%yorigin, & + DZERO, & + sin(dis%angrot * DPIO180), & + cos(dis%angrot * DPIO180), & + invert=.true.) + x = xout + y = yout + end if + ! Set status if (particle%istatus .lt. 0) then status = 1 diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index fe69f888cec..c04668dd160 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -2,7 +2,7 @@ module PrtPrpModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: DZERO, DEM1, DONE, LENFTYPE, LINELENGTH, & LENBOUNDNAME, LENPAKLOC, TABLEFT, TABCENTER, & - MNORMAL, DSAME, DEP3, DEP9 + MNORMAL, DSAME, DEP3, DEP9, DPIO180 use BndModule, only: BndType use ObsModule, only: DefaultObsIdProcessor use TableModule, only: TableType, table_cr @@ -18,10 +18,11 @@ module PrtPrpModule store_warning use SimVariablesModule, only: errmsg, warnmsg use TrackControlModule, only: TrackControlType - use GeomUtilModule, only: point_in_polygon, get_ijk, get_jk + use GeomUtilModule, only: point_in_polygon, get_ijk, get_jk, transform use MemoryManagerModule, only: mem_allocate, mem_deallocate, & mem_reallocate use ReleaseScheduleModule, only: ReleaseScheduleType, create_release_schedule + use BaseDisModule, only: DisBaseType use DisModule, only: DisType use DisvModule, only: DisvType use ErrorUtilModule, only: pstop @@ -48,6 +49,7 @@ module PrtPrpModule integer(I4B), pointer :: nparticles => null() !< number of particles released integer(I4B), pointer :: istopweaksink => null() !< weak sink option: 0 = no stop, 1 = stop integer(I4B), pointer :: istopzone => null() !< optional stop zone number: 0 = no stop zone + integer(I4B), pointer :: icoords => null() !< coordinate system option: 0 = model, 1 = global, 2 = local (unit-normalized, only for dis grids), 3 = local offset integer(I4B), pointer :: idrape => null() !< drape option: 0 = do not drape, 1 = drape to topmost active cell integer(I4B), pointer :: itrkout => null() !< binary track file integer(I4B), pointer :: itrkhdr => null() !< track header file @@ -157,6 +159,7 @@ subroutine prp_da(this) call mem_deallocate(this%stoptraveltime) call mem_deallocate(this%istopweaksink) call mem_deallocate(this%istopzone) + call mem_deallocate(this%icoords) call mem_deallocate(this%idrape) call mem_deallocate(this%nreleasepoints) call mem_deallocate(this%nreleasetimes) @@ -246,6 +249,7 @@ subroutine prp_allocate_scalars(this) call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath) call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath) call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath) + call mem_allocate(this%icoords, 'ICOORDS', this%memoryPath) call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath) call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath) call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath) @@ -269,6 +273,7 @@ subroutine prp_allocate_scalars(this) this%stoptraveltime = huge(1d0) this%istopweaksink = 0 this%istopzone = 0 + this%icoords = 0 this%idrape = 0 this%nreleasepoints = 0 this%nreleasetimes = 0 @@ -415,9 +420,12 @@ subroutine release(this, ip, trelease) ! local integer(I4B) :: irow, icol, ilay, icpl integer(I4B) :: ic, icu + integer(I4B) :: i, j, k + integer(I4B) :: icoords integer(I4B) :: np - real(DP) :: x, y, z + real(DP) :: x, y, z, xout, yout, zout real(DP) :: top, bot, hds + real(DP), allocatable :: polyverts(:, :) type(ParticleType), pointer :: particle ! Increment cumulative particle count @@ -428,9 +436,52 @@ subroutine release(this, ip, trelease) ic = this%rptnode(ip) icu = this%dis%get_nodeuser(ic) - ! Load coordinates and transform if needed + ! Load x/y coordinates and transform if needed + icoords = 0 x = this%rptx(ip) y = this%rpty(ip) + if (this%icoords == 1 .and. ( & + this%dis%angrot /= DZERO .or. & + this%dis%xorigin /= DZERO .or. & + this%dis%yorigin /= DZERO)) then + icoords = 1 + ! Transform x and y from global to model coordinates + call transform(x, y, z, & + xout, yout, zout, & + this%dis%xorigin, & + this%dis%yorigin, & + DZERO, & + sin(this%dis%angrot * DPIO180), & + cos(this%dis%angrot * DPIO180)) + x = xout + y = yout + else if (this%icoords == 2) then + ! Transform x and y from local to model coordinates + select type (dis => this%fmi%dis) + type is (DisType) + call dis%get_polyverts(ic, polyverts) + call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, i, j, k) + x = minval(polyverts(1, :)) + x * dis%delc(j) + y = minval(polyverts(2, :)) + y * dis%delr(i) + type is (DisvType) + errmsg = 'Error: LOCAL_XY may only be used on structured grids.' + call store_error(errmsg, terminate=.false.) + call store_error_unit(this%inunit, terminate=.true.) + end select + else if (this%icoords == 3) then + ! Transform x and y from local offsets to model coordinates + select type (dis => this%fmi%dis) + type is (DisType) + call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, i, j, k) + x = dis%cellx(j) + x + y = dis%celly(i) + y + type is (DisvType) + x = dis%cellxy(ic, 1) + x + y = dis%cellxy(ic, 2) + y + end select + end if + + ! Load z coordinate and transform if needed if (this%ilocalz > 0) then top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) @@ -494,6 +545,7 @@ subroutine release(this, ip, trelease) particle%ifrctrn = this%ifrctrn particle%iexmeth = this%iexmeth particle%iextend = this%iextend + particle%icoords = icoords particle%extol = this%extol ! Store the particle's state and deallocate it @@ -732,6 +784,12 @@ subroutine prp_options(this, option, found) case ('LOCAL_Z') this%ilocalz = 1 found = .true. + case ('LOCAL_XY') + this%icoords = 2 + found = .true. + case ('LOCAL_XY_OFFSET') + this%icoords = 3 + found = .true. case ('EXTEND_TRACKING') this%iextend = 1 found = .true. @@ -764,6 +822,10 @@ subroutine prp_options(this, option, found) call store_error('DEV_EXIT_SOLVE_METHOD MUST BE & &1 (BRENT) OR 2 (CHANDRUPATLA)') found = .true. + case ('DEV_GLOBAL_XY') + call this%parser%DevOpt() + this%icoords = 1 + found = .true. case default found = .false. end select @@ -909,8 +971,19 @@ subroutine prp_read_packagedata(this) y(n) = this%parser%GetDouble() z(n) = this%parser%GetDouble() + ! check coordinates + ! TODO: factor out a routine? + if (this%icoords == 2) then + if (x(n) < 0 .or. x(n) > 1) then + call store_error('Local x coordinate must fall in the unit interval') + cycle + else if (x(n) < 0 .or. x(n) > 1) then + call store_error('Local x coordinate must fall in the unit interval') + cycle + end if + end if if (this%ilocalz > 0 .and. (z(n) < 0 .or. z(n) > 1)) then - call store_error('Local z coordinate must fall in the interval [0, 1]') + call store_error('Local z coordinate must fall in the unit interval') cycle end if diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index 38049b5d51a..47957873651 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -180,8 +180,9 @@ subroutine save(this, particle, reason) end if ! Save the particle's state to any registered tracking output files - call this%trackctl%save(particle, kper=per, & - kstp=stp, reason=reason) + call this%trackctl%save(particle, this%fmi%dis, & + kper=per, kstp=stp, & + reason=reason) end subroutine save !> @brief Update particle state and check termination conditions diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index 71fc327fef7..ccd2c5a9949 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -63,6 +63,7 @@ subroutine load_mcp(this, particle, next_level, submethod) call this%load_subcell(particle, subcell) end select call method_subcell_plck%init( & + fmi=this%fmi, & cell=this%cell, & subcell=this%subcell, & trackctl=this%trackctl, & diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index 8fe7e41525f..439fba6b34d 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -216,7 +216,7 @@ subroutine apply_mcpq(this, particle, tmax) ! pass it vertically and instantaneously to the top if (particle%z > cell%defn%top) then particle%z = cell%defn%top - call this%save(particle, reason=1) ! reason=1: cell transition + call this%save(particle, reason=1) end if ! Transform particle location into local cell coordinates diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index fac00b00845..5799f96d782 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -172,7 +172,6 @@ subroutine apply_mct(this, particle, tmax) type is (CellPolyType) ! Number of vertices this%nverts = cell%defn%npolyverts - ! (Re)allocate type-bound arrays if (allocated(this%xvert)) then deallocate (this%xvert) deallocate (this%yvert) diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index 3bb1471345c..1dba3dff82b 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -20,9 +20,9 @@ module ParticleModule !! is solved for each time step. !! !! Particle coordinates may be local to the cell or - !! global/model. Routines are provided to convert a - !! particle's global coordinates to/from cell-local - !! coordinates for tracking through cell subdomains. + !! model coords. Routines are provided to convert a + !! particle's model coordinates to/from cell-local + !! coordinates for tracking through subdomains. !! !! Particles are identified by composite key, i.e., !! combinations of properties imdl, iprp, irpt, and @@ -64,6 +64,7 @@ module ParticleModule integer(I4B), public :: ifrctrn !< whether to force solving the particle with the ternary method integer(I4B), public :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method integer(I4B), public :: iextend !< whether to extend tracking beyond the end of the simulation + integer(I4B), public :: icoords !< coordinate system to use: 0 = model, 1: global contains procedure, public :: get_model_coords procedure, public :: load_particle @@ -73,31 +74,32 @@ module ParticleModule !> @brief Structure of arrays to store particles. type ParticleStoreType private - ! identity + ! identity (composite key fields) character(len=LENBOUNDNAME), dimension(:), pointer, public, contiguous :: name !< optional particle label integer(I4B), dimension(:), pointer, public, contiguous :: imdl !< index of model particle originated in integer(I4B), dimension(:), pointer, public, contiguous :: iprp !< index of release package the particle originated in integer(I4B), dimension(:), pointer, public, contiguous :: irpt !< index of release point in the particle release package the particle originated in - ! stopping criteria + ! options + integer(I4B), dimension(:), pointer, public, contiguous :: icoords !< coordinate system to use (0: model, 1: global) + integer(I4B), dimension(:), pointer, public, contiguous :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method + integer(I4B), dimension(:), pointer, public, contiguous :: iextend !< whether to extend tracking beyond the end of the simulation + integer(I4B), dimension(:), pointer, public, contiguous :: ifrctrn !< force ternary method integer(I4B), dimension(:), pointer, public, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number + real(DP), dimension(:), pointer, public, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method ! state - integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy - integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries + integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< domain hierarchy indices + integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< domain boundary indices integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user, not reduced) - integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer - integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number - integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status - real(DP), dimension(:), pointer, public, contiguous :: x !< model x coord of particle - real(DP), dimension(:), pointer, public, contiguous :: y !< model y coord of particle - real(DP), dimension(:), pointer, public, contiguous :: z !< model z coord of particle - real(DP), dimension(:), pointer, public, contiguous :: trelease !< particle release time - real(DP), dimension(:), pointer, public, contiguous :: tstop !< particle stop time - real(DP), dimension(:), pointer, public, contiguous :: ttrack !< current tracking time - integer(I4B), dimension(:), pointer, public, contiguous :: ifrctrn !< force ternary method - integer(I4B), dimension(:), pointer, public, contiguous :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method - real(DP), dimension(:), pointer, public, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method - integer(LGP), dimension(:), pointer, public, contiguous :: extend !< whether to extend tracking beyond the end of the simulation + integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer number + integer(I4B), dimension(:), pointer, public, contiguous :: izone !< zone number + integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< tracking status + real(DP), dimension(:), pointer, public, contiguous :: x !< x coordinate (model system) + real(DP), dimension(:), pointer, public, contiguous :: y !< y coordinate (model system) + real(DP), dimension(:), pointer, public, contiguous :: z !< z coordinate (model system) + real(DP), dimension(:), pointer, public, contiguous :: trelease !< release time + real(DP), dimension(:), pointer, public, contiguous :: tstop !< stop time + real(DP), dimension(:), pointer, public, contiguous :: ttrack !< tracking time contains procedure, public :: deallocate procedure, public :: num_stored @@ -141,7 +143,8 @@ subroutine allocate_particle_store(this, np, mempath) call mem_allocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_allocate(this%iexmeth, np, 'PLIEXMETH', mempath) call mem_allocate(this%extol, np, 'PLEXTOL', mempath) - call mem_allocate(this%extend, np, 'PLIEXTEND', mempath) + call mem_allocate(this%iextend, np, 'PLIEXTEND', mempath) + call mem_allocate(this%icoords, np, 'PLICOORDS', mempath) call mem_allocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath) call mem_allocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath) end subroutine allocate_particle_store @@ -170,7 +173,8 @@ subroutine deallocate (this, mempath) call mem_deallocate(this%ifrctrn, 'PLIFRCTRN', mempath) call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath) call mem_deallocate(this%extol, 'PLEXTOL', mempath) - call mem_deallocate(this%extend, 'PLIEXTEND', mempath) + call mem_deallocate(this%iextend, 'PLIEXTEND', mempath) + call mem_deallocate(this%icoords, 'PLICOORDS', mempath) call mem_deallocate(this%idomain, 'PLIDOMAIN', mempath) call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath) end subroutine deallocate @@ -202,7 +206,8 @@ subroutine resize(this, np, mempath) call mem_reallocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath) call mem_reallocate(this%extol, np, 'PLEXTOL', mempath) - call mem_reallocate(this%extend, np, 'PLIEXTEND', mempath) + call mem_reallocate(this%iextend, np, 'PLIEXTEND', mempath) + call mem_reallocate(this%icoords, np, 'PLICOORDS', mempath) call mem_reallocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath) call mem_reallocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath) end subroutine resize @@ -246,7 +251,8 @@ subroutine load_particle(this, store, imdl, iprp, ip) this%ifrctrn = store%ifrctrn(ip) this%iexmeth = store%iexmeth(ip) this%extol = store%extol(ip) - this%iextend = store%extend(ip) + this%iextend = store%iextend(ip) + this%icoords = store%icoords(ip) end subroutine load_particle !> @brief Save a particle's state to the particle store @@ -282,7 +288,8 @@ subroutine save_particle(this, particle, ip) this%ifrctrn(ip) = particle%ifrctrn this%iexmeth(ip) = particle%iexmeth this%extol(ip) = particle%extol - this%extend(ip) = particle%iextend + this%iextend(ip) = particle%iextend + this%icoords(ip) = particle%icoords end subroutine save_particle !> @brief Apply the given global-to-local transformation to the particle. @@ -306,7 +313,6 @@ subroutine transform_coords(this, xorigin, yorigin, zorigin, & this%zorigin = DZERO this%sinrot = DZERO this%cosrot = DONE - this%cosrot = DONE this%transformed = .false. return end if @@ -333,7 +339,7 @@ subroutine transform_coords(this, xorigin, yorigin, zorigin, & this%transformed = .true. end subroutine transform_coords - !> @brief Return the particle's model (global) coordinates. + !> @brief Return the particle's model coordinates. subroutine get_model_coords(this, x, y, z) use GeomUtilModule, only: transform, compose class(ParticleType), intent(inout) :: this !< particle