From 5ce7eadb724adaefbbff4ccd228ea367e22a3bb0 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 14 Dec 2023 15:13:38 -0500 Subject: [PATCH] more cleanup, remove unused code --- src/Model/ParticleTracking/prt1.f90 | 2 +- src/Model/ParticleTracking/prt1prp1.f90 | 16 +++--- src/Solution/ParticleTracker/Method.f90 | 8 +-- .../ParticleTracker/MethodCellPassToBot.f90 | 2 +- .../ParticleTracker/MethodCellPollock.f90 | 13 +++-- .../ParticleTracker/MethodCellPollockQuad.f90 | 51 +++++++++---------- .../ParticleTracker/MethodCellTernary.f90 | 33 ++++++------ src/Solution/ParticleTracker/MethodDis.f90 | 18 +++---- src/Solution/ParticleTracker/MethodDisv.f90 | 22 ++++---- .../ParticleTracker/MethodSubcellPollock.f90 | 6 +-- .../ParticleTracker/MethodSubcellTernary.f90 | 8 +-- src/Solution/ParticleTracker/Particle.f90 | 50 +++++++++--------- .../ParticleTracker/TernarySolveTrack.f90 | 51 ------------------- 13 files changed, 112 insertions(+), 168 deletions(-) diff --git a/src/Model/ParticleTracking/prt1.f90 b/src/Model/ParticleTracking/prt1.f90 index eef14fba1e8..c0c2e72247a 100644 --- a/src/Model/ParticleTracking/prt1.f90 +++ b/src/Model/ParticleTracking/prt1.f90 @@ -438,7 +438,7 @@ subroutine prt_cq_sto(this) ! refine these conditions as necessary ! (status 8 is permanently unreleased) if ((istatus > 0) .and. (istatus /= 8)) then - n = packobj%particles%iTrackingDomain(np, 2) + n = packobj%particles%idomain(np, 2) ! -- Each particle currently assigned unit mass this%masssto(n) = this%masssto(n) + DONE end if diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index b335c2599f3..a098f37ee27 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -500,14 +500,14 @@ subroutine prp_ad(this) particle%trelease = trelease particle%tstop = tstop particle%ttrack = trelease - particle%iTrackingDomain(0) = 0 - particle%iTrackingDomainBoundary(0) = 0 - particle%iTrackingDomain(1) = 0 - particle%iTrackingDomainBoundary(1) = 0 - particle%iTrackingDomain(2) = ic - particle%iTrackingDomainBoundary(2) = 0 - particle%iTrackingDomain(3) = 0 - particle%iTrackingDomainBoundary(3) = 0 + particle%idomain(0) = 0 + particle%iboundary(0) = 0 + particle%idomain(1) = 0 + particle%iboundary(1) = 0 + particle%idomain(2) = ic + particle%iboundary(2) = 0 + particle%idomain(3) = 0 + particle%iboundary(3) = 0 ! -- Add particle to particle list call this%particles%update_from_particle(particle, np) diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index ac83f5ef6e6..1ea7c94de34 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -112,22 +112,22 @@ recursive subroutine track(this, particle, level, tmax) end subroutine track !> @brief Try passing the particle to the next subdomain - subroutine try_pass(this, particle, levelNext, advancing) + subroutine try_pass(this, particle, next_level, advancing) class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle - integer :: levelNext + integer :: next_level logical :: advancing ! tracking submethod marked tracking complete? ! reset domain boundary flag and don't advance if (.not. particle%advancing) then - particle%iTrackingDomainBoundary = 0 + particle%iboundary = 0 advancing = .false. else ! otherwise pass particle to next subdomain ! and if it's on a boundary, stop advancing call this%pass(particle) - if (particle%iTrackingDomainBoundary(levelNext - 1) .ne. 0) & + if (particle%iboundary(next_level - 1) .ne. 0) & advancing = .false. end if end subroutine try_pass diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 index 9863e90be2d..4014ad45eab 100644 --- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 +++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 @@ -53,7 +53,7 @@ subroutine apply_ptb(this, particle, tmax) call this%update(particle, this%defn) if (.not. particle%advancing) return particle%z = this%defn%bot - particle%iTrackingDomainBoundary(2) = this%defn%npolyverts + 2 + particle%iboundary(2) = this%defn%npolyverts + 2 call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=1) ! reason=1: cell transition end subroutine apply_ptb diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index c79b2f72929..b1851c6da44 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -59,10 +59,11 @@ subroutine load_mcp(this, particle, next_level, submethod) select type (subcell => this%subcell) type is (SubcellRectType) - call this%load_subcell(particle, next_level, subcell) + call this%load_subcell(particle, subcell) end select call method_subcell_plck%init(subcell=this%subcell, trackctl=this%trackctl) submethod => method_subcell_plck + particle%idomain(next_level) = 1 end subroutine load_mcp !> @brief Having exited the lone subcell, pass the particle to the cell face @@ -74,7 +75,7 @@ subroutine pass_mcp(this, particle) ! -- local integer :: exitFace, inface - exitFace = particle%iTrackingDomainBoundary(3) + exitFace = particle%iboundary(3) ! -- Map subcell exit face to cell face select case (exitFace) ! note: exitFace uses Dave's iface convention case (0) @@ -93,7 +94,7 @@ subroutine pass_mcp(this, particle) inface = 7 end select if (inface .eq. -1) then - particle%iTrackingDomainBoundary(2) = 0 + particle%iboundary(2) = 0 else if ((inface .ge. 1) .and. (inface .le. 4)) then ! -- Account for local cell rotation @@ -103,7 +104,7 @@ subroutine pass_mcp(this, particle) end select if (inface .gt. 4) inface = inface - 4 end if - particle%iTrackingDomainBoundary(2) = inface + particle%iboundary(2) = inface end if end subroutine pass_mcp @@ -160,18 +161,16 @@ end subroutine apply_mcp !> @brief Loads the lone rectangular subcell from the rectangular cell !! kluge note: is levelNext needed here and in similar "load" routines??? - subroutine load_subcell(this, particle, next_level, subcell) ! + subroutine load_subcell(this, particle, subcell) ! ! -- dummy class(MethodCellPollockType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle - integer, intent(in) :: next_level type(SubcellRectType), intent(inout) :: subcell select type (cell => this%cell) type is (CellRectType) ! -- Set subcell number to 1 subcell%isubcell = 1 - particle%iTrackingDomain(next_level) = 1 ! kluge note: is this the place to set this??? ! -- Subcell calculations will be done in local subcell coordinates subcell%dx = cell%dx diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index 19e44f504d4..fc80a68f9ea 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -59,7 +59,7 @@ subroutine load_mcpq(this, particle, next_level, submethod) select type (subcell => this%subcell) type is (SubcellRectType) - call this%load_subcell(particle, next_level, subcell) + call this%load_subcell(particle, subcell) end select call method_subcell_plck%init(subcell=this%subcell, trackctl=this%trackctl) submethod => method_subcell_plck @@ -75,8 +75,8 @@ subroutine pass_mcpq(this, particle) select type (cell => this%cell) type is (CellRectQuadType) - exitFace = particle%iTrackingDomainBoundary(3) - isc = particle%iTrackingDomain(3) + exitFace = particle%iboundary(3) + isc = particle%idomain(3) npolyverts = cell%defn%npolyverts select case (exitFace) ! kluge note: exitFace uses Dave's iface convention @@ -87,13 +87,13 @@ subroutine pass_mcpq(this, particle) select case (isc) case (1) ! -- W face, subcell 1 --> E face, subcell 4 (cell interior) - particle%iTrackingDomain(3) = 4 - particle%iTrackingDomainBoundary(3) = 2 + particle%idomain(3) = 4 + particle%iboundary(3) = 2 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (2) ! -- W face, subcell 2 --> E face, subcell 3 (cell interior) - particle%iTrackingDomain(3) = 3 - particle%iTrackingDomainBoundary(3) = 2 + particle%idomain(3) = 3 + particle%iboundary(3) = 2 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (3) ! -- W face, subcell 3 (cell face) @@ -116,21 +116,21 @@ subroutine pass_mcpq(this, particle) infaceoff = -1 case (3) ! -- E face, subcell 3 --> W face, subcell 2 (cell interior) - particle%iTrackingDomain(3) = 2 - particle%iTrackingDomainBoundary(3) = 1 + particle%idomain(3) = 2 + particle%iboundary(3) = 1 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (4) ! -- E face, subcell 4 --> W face subcell 1 (cell interior) - particle%iTrackingDomain(3) = 1 - particle%iTrackingDomainBoundary(3) = 1 + particle%idomain(3) = 1 + particle%iboundary(3) = 1 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 end select case (3) select case (isc) case (1) ! -- S face, subcell 1 --> N face, subcell 2 (cell interior) - particle%iTrackingDomain(3) = 2 - particle%iTrackingDomainBoundary(3) = 4 + particle%idomain(3) = 2 + particle%iboundary(3) = 4 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (2) ! -- S face, subcell 2 (cell face) @@ -142,8 +142,8 @@ subroutine pass_mcpq(this, particle) infaceoff = -1 case (4) ! -- S face, subcell 4 --> N face, subcell 3 (cell interior) - particle%iTrackingDomain(3) = 3 - particle%iTrackingDomainBoundary(3) = 4 + particle%idomain(3) = 3 + particle%iboundary(3) = 4 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 end select case (4) @@ -154,13 +154,13 @@ subroutine pass_mcpq(this, particle) infaceoff = -1 case (2) ! -- N face, subcell 2 --> S face, subcell 1 (cell interior) - particle%iTrackingDomain(3) = 1 - particle%iTrackingDomainBoundary(3) = 3 + particle%idomain(3) = 1 + particle%iboundary(3) = 3 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (3) ! -- N face, subcell 3 --> S face, subcell 4 (cell interior) - particle%iTrackingDomain(3) = 4 - particle%iTrackingDomainBoundary(3) = 3 + particle%idomain(3) = 4 + particle%iboundary(3) = 3 inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0 case (4) ! -- N face, subcell 4 (cell face) @@ -175,9 +175,9 @@ subroutine pass_mcpq(this, particle) inface = npolyverts + 3 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface end select if (inface .eq. -1) then - particle%iTrackingDomainBoundary(2) = 0 + particle%iboundary(2) = 0 else if (inface .eq. 0) then - particle%iTrackingDomainBoundary(2) = 0 + particle%iboundary(2) = 0 else if ((inface .ge. 1) .and. (inface .le. 4)) then ! -- Account for local cell rotation @@ -186,7 +186,7 @@ subroutine pass_mcpq(this, particle) inface = cell%irectvert(inface) + infaceoff if (inface .lt. 1) inface = inface + npolyverts end if - particle%iTrackingDomainBoundary(2) = inface + particle%iboundary(2) = inface end if end select end subroutine pass_mcpq @@ -236,11 +236,10 @@ subroutine apply_mcpq(this, particle, tmax) end subroutine apply_mcpq !> @brief Load the rectangular subcell from the rectangular cell - subroutine load_subcell(this, particle, next_level, subcell) + subroutine load_subcell(this, particle, subcell) ! -- dummy class(MethodCellPollockQuadType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle - integer, intent(in) :: next_level class(SubcellRectType), intent(inout) :: subcell ! -- local double precision :: dx, dy, dz, areax, areay, areaz @@ -255,7 +254,7 @@ subroutine load_subcell(this, particle, next_level, subcell) factor = factor / cell%defn%porosity npolyverts = cell%defn%npolyverts - isc = particle%iTrackingDomain(3) + isc = particle%idomain(3) ! -- Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D, ! -- respectively @@ -290,7 +289,7 @@ subroutine load_subcell(this, particle, next_level, subcell) subcell%isubcell = isc ! kluge note: as a matter of form, do we want to allow ! this subroutine to modify the particle??? - particle%iTrackingDomain(3) = isc + particle%idomain(3) = isc ! kluge note: initial insubface is not currently being determined end if dx = 5d-1 * dx diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index 32a48b3fd6a..7b1c8800853 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -67,7 +67,7 @@ subroutine load_mct(this, particle, next_level, submethod) select type (subcell => this%subcell) type is (SubcellTriType) - call this%load_subcell(particle, next_level, subcell) + call this%load_subcell(particle, subcell) end select call method_subcell_tern%init(subcell=this%subcell, trackctl=this%trackctl) submethod => method_subcell_tern @@ -81,8 +81,8 @@ subroutine pass_mct(this, particle) ! local integer :: isc, exitFace, inface, npolyverts - exitFace = particle%iTrackingDomainBoundary(3) - isc = particle%iTrackingDomain(3) + exitFace = particle%iboundary(3) + isc = particle%idomain(3) select type (cell => this%cell) type is (CellPolyType) npolyverts = cell%defn%npolyverts @@ -100,15 +100,15 @@ subroutine pass_mct(this, particle) ! -- Subcell face --> next subcell in "cycle" (cell interior) isc = isc + 1 if (isc .gt. npolyverts) isc = 1 - particle%iTrackingDomain(3) = isc - particle%iTrackingDomainBoundary(3) = 3 + particle%idomain(3) = isc + particle%iboundary(3) = 3 inface = 0 case (3) ! -- Subcell face --> preceding subcell in "cycle" (cell interior) isc = isc - 1 if (isc .lt. 1) isc = npolyverts - particle%iTrackingDomain(3) = isc - particle%iTrackingDomainBoundary(3) = 2 + particle%idomain(3) = isc + particle%iboundary(3) = 2 inface = 0 case (4) ! -- Subcell bottom (cell bottom) @@ -118,11 +118,11 @@ subroutine pass_mct(this, particle) inface = npolyverts + 3 end select if (inface .eq. -1) then - particle%iTrackingDomainBoundary(2) = 0 + particle%iboundary(2) = 0 else if (inface .eq. 0) then - particle%iTrackingDomainBoundary(2) = 0 + particle%iboundary(2) = 0 else - particle%iTrackingDomainBoundary(2) = inface + particle%iboundary(2) = inface end if end subroutine pass_mct @@ -228,14 +228,13 @@ subroutine apply_mct(this, particle, tmax) end select end subroutine apply_mct - !> @brief Loads the triangular subcell from the polygonal cell - subroutine load_subcell(this, particle, next_level, subcell) + !> @brief Loads a triangular subcell from the polygonal cell + subroutine load_subcell(this, particle, subcell) ! -- modules use ParticleModule, only: get_particle_id ! -- dummy class(MethodCellTernaryType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle - integer, intent(in) :: next_level class(SubcellTriType), intent(inout) :: subcell ! -- local integer :: ic, isc, npolyverts @@ -249,10 +248,8 @@ subroutine load_subcell(this, particle, next_level, subcell) type is (CellPolyType) ic = cell%defn%icell subcell%icell = ic - isc = particle%iTrackingDomain(3) + isc = particle%idomain(3) npolyverts = cell%defn%npolyverts - - ! -- Find subcell if not known ! kluge note: from "find_init_triangle", todo: move there if (isc .le. 0) then xi = particle%x yi = particle%y @@ -279,7 +276,7 @@ subroutine load_subcell(this, particle, next_level, subcell) d01 = x0 * y1rel - y0 * x1rel alphai = (di2 - d02) / d12 betai = -(di1 - d01) / d12 - ! kluge note: can iTrackingDomainBoundary(2) be used to identify the subcell? + ! kluge note: can iboundary(2) be used to identify the subcell? betatol = -1e-7 ! kluge ! kluge note: think this handles points on triangle boundaries ok if ((alphai .ge. 0d0) .and. & @@ -297,7 +294,7 @@ subroutine load_subcell(this, particle, next_level, subcell) ! subcellTri%isubcell = isc ! kluge note: as a matter of form, do we want to allow ! this subroutine to modify the particle??? - particle%iTrackingDomain(3) = isc + particle%idomain(3) = isc end if end if subcell%isubcell = isc diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index 05add8ce162..fde9e82057c 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -73,7 +73,7 @@ subroutine load_dis(this, particle, next_level, submethod) type is (CellRectType) select type (dis => this%fmi%dis) type is (GwfDisType) - ic = particle%iTrackingDomain(next_level) + ic = particle%idomain(next_level) call this%load_cell_defn(ic, cell%defn) ! -- If cell is active but dry, select and initialize @@ -132,7 +132,7 @@ subroutine pass_dis(this, particle) integer :: inface, ipos, ic, icu, inbr, idiag, ilay, irow, icol real(DP) :: z, zrel, topfrom, botfrom, top, bot, sat - inface = particle%iTrackingDomainBoundary(2) + inface = particle%iboundary(2) z = particle%z select type (cell => this%cell) @@ -142,20 +142,20 @@ subroutine pass_dis(this, particle) inbr = cell%defn%facenbr(inface) if (inbr .eq. 0) then ! -- Exterior face; no neighbor to map to - ! particle%iTrackingDomain(1) = 0 - ! particle%iTrackingDomain(2) = 0 ! kluge note: set a "has_exited" attribute instead??? - ! particle%iTrackingDomain(1) = -abs(particle%iTrackingDomain(1)) ! kluge??? - ! particle%iTrackingDomain(2) = -abs(particle%iTrackingDomain(2)) ! kluge??? + ! particle%idomain(1) = 0 + ! particle%idomain(2) = 0 ! kluge note: set a "has_exited" attribute instead??? + ! particle%idomain(1) = -abs(particle%idomain(1)) ! kluge??? + ! particle%idomain(2) = -abs(particle%idomain(2)) ! kluge??? particle%istatus = 2 ! kluge note: use -2 to allow check for transfer to another model??? particle%advancing = .false. call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=3) ! reason=3: termination - ! particle%iTrackingDomainBoundary(2) = -1 + ! particle%iboundary(2) = -1 else idiag = dis%con%ia(cell%defn%icell) ipos = idiag + inbr ic = dis%con%ja(ipos) ! kluge note: use PRT model's DIS instead of fmi's??? - particle%iTrackingDomain(2) = ic + particle%idomain(2) = ic ! compute and set user node number and layer on particle icu = dis%get_nodeuser(ic) @@ -178,7 +178,7 @@ subroutine pass_dis(this, particle) else if (inface .eq. 7) then inface = 6 end if - particle%iTrackingDomainBoundary(2) = inface + particle%iboundary(2) = inface if (inface < 5) then ! -- Map z between cells topfrom = cell%defn%top diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index ec5b76e462a..a9444a109aa 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -80,7 +80,7 @@ subroutine load_disv(this, particle, next_level, submethod) select type (cell => this%cell) type is (CellPolyType) ! load cell definition - ic = particle%iTrackingDomain(next_level) ! kluge note: is cell number always known coming in? + ic = particle%idomain(next_level) ! kluge note: is cell number always known coming in? call this%load_cell_defn(ic, cell%defn) 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 @@ -119,7 +119,7 @@ subroutine pass_disv(this, particle) integer :: inface, ipos, ic, icu, inbr, idiag, icpl, ilay double precision :: z - inface = particle%iTrackingDomainBoundary(2) + inface = particle%iboundary(2) z = particle%z select type (cell => this%cell) @@ -129,20 +129,20 @@ subroutine pass_disv(this, particle) inbr = cell%defn%facenbr(inface) if (inbr .eq. 0) then ! -- Exterior face; no neighbor to map to - ! particle%iTrackingDomain(1) = 0 - ! particle%iTrackingDomain(2) = 0 ! kluge note: "has_exited" attribute instead??? - ! particle%iTrackingDomain(1) = -abs(particle%iTrackingDomain(1)) ! kluge??? - ! particle%iTrackingDomain(2) = -abs(particle%iTrackingDomain(2)) ! kluge??? + ! particle%idomain(1) = 0 + ! particle%idomain(2) = 0 ! kluge note: "has_exited" attribute instead??? + ! particle%idomain(1) = -abs(particle%idomain(1)) ! kluge??? + ! particle%idomain(2) = -abs(particle%idomain(2)) ! kluge??? particle%istatus = 2 ! kluge note, todo: use -2 to check for transfer to another model??? particle%advancing = .false. call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=3) ! reason=3: termination - ! particle%iTrackingDomainBoundary(2) = -1 + ! particle%iboundary(2) = -1 else idiag = dis%con%ia(cell%defn%icell) ipos = idiag + inbr ic = dis%con%ja(ipos) ! kluge note, todo: use PRT model's DIS instead of fmi's?? - particle%iTrackingDomain(2) = ic + particle%idomain(2) = ic ! compute and set user node number and layer on particle icu = dis%get_nodeuser(ic) @@ -151,9 +151,9 @@ subroutine pass_disv(this, particle) particle%ilay = ilay call this%map_neighbor(cell%defn, inface, z) - particle%iTrackingDomainBoundary(2) = inface - particle%iTrackingDomain(3:) = 0 - particle%iTrackingDomainBoundary(3:) = 0 + particle%iboundary(2) = inface + particle%idomain(3:) = 0 + particle%iboundary(3:) = 0 particle%z = z ! -- Update cell-cell flows of particle mass. ! Every particle is currently assigned unit mass. diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index f53274beadc..a7a3f6d2c20 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -125,14 +125,14 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- Check for no exit face if ((statusVX .eq. 3) .and. (statusVY .eq. 3) .and. (statusVZ .eq. 3)) then ! exitFace = 0 - ! particle%iTrackingDomainBoundary(3) = exitFace + ! 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%iTrackingDomain(2) + "cell", particle%idomain(2) call pstop(1) end if @@ -218,7 +218,7 @@ subroutine track_subcell(this, subcell, particle, tmax) particle%y = y * subcell%dy particle%z = z * subcell%dz particle%ttrack = t - particle%iTrackingDomainBoundary(3) = exitFace + particle%iboundary(3) = exitFace ! -- Save particle track record if (reason > -1) & diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 8141f998914..b449d1397c6 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -138,7 +138,7 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- Traverse triangular subcell ntdebug = -999 ! kluge debug bludebug - itrifaceenter = particle%iTrackingDomainBoundary(3) - 1 + itrifaceenter = particle%iboundary(3) - 1 if (itrifaceenter .eq. -1) itrifaceenter = 999 ! kluge note: can probably avoid calculating alpexit @@ -156,13 +156,13 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- Check for no exit face if ((itopbotexit .eq. 0) .and. (itrifaceexit .eq. 0)) then ! exitFace = 0 - ! particle%iTrackingDomainBoundary(3) = exitFace + ! 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%iTrackingDomain(2) + "cell", particle%idomain(2) call pstop(1) end if @@ -252,7 +252,7 @@ subroutine track_subcell(this, subcell, particle, tmax) particle%y = y particle%z = z particle%ttrack = t - particle%iTrackingDomainBoundary(3) = exitFace + particle%iboundary(3) = exitFace ! -- Save particle track record if (reason > -1) & diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index d1c312b1640..97878ec6262 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -37,8 +37,8 @@ module ParticleModule integer(I4B), public :: istopzone !< stop zone number ! state - integer(I4B), allocatable, public :: iTrackingDomain(:) !< array of indices for domains in the tracking domain hierarchy - integer(I4B), allocatable, public :: iTrackingDomainBoundary(:) !< array of indices for tracking domain boundaries + integer(I4B), allocatable, public :: idomain(:) !< array of indices for domains in the tracking domain hierarchy + integer(I4B), allocatable, public :: iboundary(:) !< array of indices for tracking domain boundaries integer(I4B), public :: icu !< user cell (node) number integer(I4B), public :: ilay !< layer integer(I4B), public :: izone !< current zone number @@ -82,8 +82,8 @@ module ParticleModule integer(I4B), dimension(:), pointer, contiguous :: istopzone !< stop zone number ! state - integer(I4B), dimension(:, :), allocatable :: iTrackingDomain !< array of indices for domains in the tracking domain hierarchy - integer(I4B), dimension(:, :), allocatable :: iTrackingDomainBoundary !< array of indices for tracking domain boundaries + integer(I4B), dimension(:, :), allocatable :: idomain !< array of indices for domains in the tracking domain hierarchy + integer(I4B), dimension(:, :), allocatable :: iboundary !< array of indices for tracking domain boundaries integer(I4B), dimension(:), pointer, contiguous :: icu !< cell number (user, not reduced) integer(I4B), dimension(:), pointer, contiguous :: ilay !< layer integer(I4B), dimension(:), pointer, contiguous :: izone !< current zone number @@ -108,15 +108,15 @@ module ParticleModule subroutine create_particle(particle) type(ParticleType), pointer :: particle !< particle allocate (particle) - allocate (particle%iTrackingDomain(levelMin:levelMax)) - allocate (particle%iTrackingDomainBoundary(levelMin:levelMax)) + allocate (particle%idomain(levelMin:levelMax)) + allocate (particle%iboundary(levelMin:levelMax)) end subroutine create_particle !> @brief Destroy a particle subroutine destroy_particle(this) class(ParticleType), intent(inout) :: this !< particle - deallocate (this%iTrackingDomain) - deallocate (this%iTrackingDomainBoundary) + deallocate (this%idomain) + deallocate (this%iboundary) end subroutine destroy_particle !> @brief Allocate particle arrays @@ -136,8 +136,8 @@ subroutine allocate_arrays(this, np, lmin, lmax, mempath) call mem_allocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath) ! -- kluge todo: update mem_allocate to allow custom range of indices? ! e.g. here we want to allocate 0-4 for trackdomain levels, not 1-5 - allocate (this%iTrackingDomain(np, lmin:lmax)) - allocate (this%iTrackingDomainBoundary(np, lmin:lmax)) + allocate (this%idomain(np, lmin:lmax)) + allocate (this%iboundary(np, lmin:lmax)) call mem_allocate(this%icu, np, 'PLICU', mempath) call mem_allocate(this%ilay, np, 'PLILAY', mempath) call mem_allocate(this%izone, np, 'PLIZONE', mempath) @@ -165,8 +165,8 @@ subroutine deallocate_arrays(this, mempath) call mem_deallocate(this%iprp, 'PLIPRP', mempath) call mem_deallocate(this%irpt, 'PLIRPT', mempath) call mem_deallocate(this%name, 'PLNAME', mempath) - deallocate (this%iTrackingDomain) - deallocate (this%iTrackingDomainBoundary) + deallocate (this%idomain) + deallocate (this%iboundary) call mem_deallocate(this%icu, 'PLICU', mempath) call mem_deallocate(this%ilay, 'PLILAY', mempath) call mem_deallocate(this%izone, 'PLIZONE', mempath) @@ -213,12 +213,12 @@ subroutine reallocate_arrays(this, np, mempath) ! resize first dimension of 2D arrays ! todo: memory manager support? call ExpandArray2D( & - this%iTrackingDomain, & - size(this%iTrackingDomain(:, 1) - np), & + this%idomain, & + size(this%idomain(:, 1) - np), & 0) call ExpandArray2D( & - this%iTrackingDomainBoundary, & - size(this%iTrackingDomainBoundary(:, 1) - np), & + this%iboundary, & + size(this%iboundary(:, 1) - np), & 0) end subroutine reallocate_arrays @@ -256,11 +256,11 @@ subroutine load_from_list(this, partlist, imdl, iprp, ip) this%tstop = partlist%tstop(ip) this%ttrack = partlist%ttrack(ip) this%advancing = .true. - this%iTrackingDomain(levelMin:levelMax) = & - partlist%iTrackingDomain(ip, levelMin:levelMax) - this%iTrackingDomain(1) = imdl - this%iTrackingDomainBoundary(levelMin:levelMax) = & - partlist%iTrackingDomainBoundary(ip, levelMin:levelMax) + this%idomain(levelMin:levelMax) = & + partlist%idomain(ip, levelMin:levelMax) + this%idomain(1) = imdl + this%iboundary(levelMin:levelMax) = & + partlist%iboundary(ip, levelMin:levelMax) end subroutine load_from_list @@ -286,14 +286,14 @@ subroutine update_from_particle(this, particle, ip) this%trelease(ip) = particle%trelease this%tstop(ip) = particle%tstop this%ttrack(ip) = particle%ttrack - this%iTrackingDomain( & + this%idomain( & ip, & levelMin:levelMax) = & - particle%iTrackingDomain(levelMin:levelMax) - this%iTrackingDomainBoundary( & + particle%idomain(levelMin:levelMax) + this%iboundary( & ip, & levelMin:levelMax) = & - particle%iTrackingDomainBoundary(levelMin:levelMax) + particle%iboundary(levelMin:levelMax) end subroutine update_from_particle diff --git a/src/Solution/ParticleTracker/TernarySolveTrack.f90 b/src/Solution/ParticleTracker/TernarySolveTrack.f90 index defc93a992c..e959862cc33 100644 --- a/src/Solution/ParticleTracker/TernarySolveTrack.f90 +++ b/src/Solution/ParticleTracker/TernarySolveTrack.f90 @@ -86,57 +86,6 @@ subroutine traverse_triangle(ntmax, nsave, diff, rdiff, & end subroutine - ! todo: reinstate - ! !> @brief Find initial triangular subcell - ! !! - ! !! kluge note: not currently used but maybe could be in load_subcell of MethodCellTernary??? - ! !! - ! !< - ! subroutine find_init_triangle(numvermx,nvert,ncpl,xi,yi,xctr,yctr,icell,itri) - ! use Ternary, only: ivert_polygon, x_vert, y_vert - ! implicit double precision(a - h, o - z) - ! ! - ! ! -- Find initial triangular subcell - ! itri = 0 - ! ! numver = num_vertices(numvermx,ncpl,icell) - ! numver = 4 ! kluge - ! do iv = 1, numver - ! iv0 = iv - ! iv1 = iv + 1 - ! if (iv1 .gt. numver) iv1 = 1 - ! ipv0 = ivert_polygon(icell, iv0) - ! ipv1 = ivert_polygon(icell, iv1) - ! x0 = x_vert(ipv0) - ! y0 = y_vert(ipv0) - ! x1 = x_vert(ipv1) - ! y1 = y_vert(ipv1) - ! x2 = xctr - ! y2 = yctr - ! x1rel = x1 - x0 - ! y1rel = y1 - y0 - ! x2rel = x2 - x0 - ! y2rel = y2 - y0 - ! di2 = xi * y2rel - yi * x2rel - ! d02 = x0 * y2rel - y0 * x2rel - ! d12 = x1rel * y2rel - y1rel * x2rel - ! di1 = xi * y1rel - yi * x1rel - ! d01 = x0 * y1rel - y0 * x1rel - ! alphai = (di2 - d02) / d12 - ! betai = -(di1 - d01) / d12 - ! if ((alphai .ge. 0d0) .and. (betai .ge. 0d0) .and. (alphai + betai .lt. 1d0)) then ! kluge note: think this handles points on triangle boundaries ok - ! itri = iv ! but maybe not!!!!!!!!!!!! - ! exit ! kluge note: doesn't handle particle smack on cell center - ! end if - ! end do - ! if (itri .eq. 0) then - ! write (*, '(A)') "error -- initial triangle not found" ! kluge - ! write (69, '(A)') "error -- initial triangle not found" - ! !!pause - ! stop - ! end if - ! - ! end subroutine - !> @brief Set coordinates to "canonical" configuration subroutine canonical(x0, y0, x1, y1, x2, y2, & v0x, v0y, v1x, v1y, v2x, v2y, &