From b8ffceeb1e1c320a905f3c0acc3243723d247744 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Mon, 23 Dec 2024 07:25:40 -0500 Subject: [PATCH] finally --- .../ParticleTracker/MethodCellTernary.f90 | 19 ++++------- .../ParticleTracker/MethodSubcellTernary.f90 | 34 ++++++++----------- .../ParticleTracker/TernarySolveTrack.f90 | 8 ----- 3 files changed, 22 insertions(+), 39 deletions(-) diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index 1385695d5ec..6b4f62a72d5 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -267,9 +267,6 @@ subroutine load_subcell(this, particle, subcell) y1 = this%yvertnext(iv0) x2 = this%xctr y2 = this%yctr - - ! print *, "subcell in load_subcell: ", ic, isc, x0, y0, x1, y1, x2, y2 - x1rel = x1 - x0 y1rel = y1 - y0 x2rel = x2 - x0 @@ -362,7 +359,7 @@ subroutine vertvelo(this) real(DP), allocatable, dimension(:) :: vm0y real(DP), allocatable, dimension(:) :: vm1x real(DP), allocatable, dimension(:) :: vm1y - real(DP), allocatable, dimension(:,:) :: poly + real(DP), allocatable, dimension(:, :) :: poly integer(I4B) :: nvert select type (cell => this%cell) @@ -413,18 +410,16 @@ subroutine vertvelo(this) this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa end if + + ! TODO: do we need to check if center is in polygon? + ! only if using the centroid method which misbehaves + ! ! allocate(poly(2, nvert)) ! poly(1,:) = this%xvert ! poly(2,:) = this%yvert ! if (.not. point_in_polygon(this%xctr, this%yctr, poly)) then - ! print *, "not in poly!!!!!! ", this%cell%defn%icell - ! if (nvert == 3) then - ! this%xctr = sum(this%xvert) / 3.0_DP - ! this%yctr = sum(this%yvert) / 3.0_DP - ! else - ! print *, "error -- centroid not in cell ", this%cell%defn%icell, this%xctr, this%yctr - ! call pstop(1) - ! end if + ! print *, "error -- centroid not in cell ", this%cell%defn%icell, this%xctr, this%yctr + ! call pstop(1) ! end if ! deallocate(poly) diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 9859961b7b5..5f368bdcff2 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -31,9 +31,9 @@ module MethodSubcellTernaryModule !> @brief Create a new ternary subcell tracking method. subroutine create_method_subcell_ternary(method) - ! -- dummy + ! dummy type(MethodSubcellTernaryType), pointer :: method - ! -- local + ! local type(SubcellTriType), pointer :: subcell allocate (method) @@ -167,8 +167,6 @@ subroutine track_subcell(this, subcell, particle, tmax) vzbot = subcell%vzbot vztop = subcell%vztop - ! print *, "subcell in track_subcell: ", subcell%icell, subcell%isubcell, x0, y0, x1, y1, x2, y2 - ! Transform coordinates to the "canonical" configuration: ! barycentric in two dimensions with alpha, beta & gamma ! such that at f2 alpha = 0, f0 beta = 0, f1 gamma = 0. @@ -223,7 +221,7 @@ subroutine track_subcell(this, subcell, particle, tmax) ! Exit through lateral subcell face exitFace = itrifaceexit dtexit = dtexitxy - else if (dtexitz < dtexitxy) then + else if (dtexitz < dtexitxy .and. dtexitz >= 0.0_DP) then ! Exit through top or bottom if (itopbotexit == -1) then exitFace = 4 @@ -231,6 +229,11 @@ subroutine track_subcell(this, subcell, particle, tmax) exitFace = 5 end if dtexit = dtexitz + else + particle%istatus = 9 + particle%advancing = .false. + call this%save(particle, reason=3) + return end if texit = particle%ttrack + dtexit t0 = particle%ttrack @@ -262,8 +265,8 @@ subroutine track_subcell(this, subcell, particle, tmax) ! (local, unscaled) and other properties. The particle may at this ! point lie on a boundary of the subcell or may still be within it. 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. + ! 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 - t0 exitFace = 0 @@ -271,8 +274,8 @@ subroutine track_subcell(this, subcell, particle, tmax) particle%advancing = .false. reason = 2 ! timestep end else - ! -- The computed exit time is less than or equal to the maximum time, - ! -- so set final time for particle trajectory equal to exit time. + ! The computed exit time is less than or equal to the maximum time, + ! so set final time for particle trajectory equal to exit time. t = texit dt = dtexit reason = 1 ! (sub)cell transition @@ -426,13 +429,6 @@ subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, & end if end if - if (vr < DZERO) then - print *, "vr v1a v2a: ", vr, v1a, v2a - print *, "vr1 vr2: ", vr1, vr2 - print *, "v1 v2: ", v1, v2 - print *, "v: ", v - end if - ! Compute travel time to exit face. Return with status = 0 dt = log(abs(vr)) / dvdx status = 0 @@ -498,13 +494,13 @@ subroutine calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, & else ! otherwise calculate z. if (izstatus .eq. 2) then - ! -- vz uniformly zero + ! vz uniformly zero z = zi else if (izstatus .eq. 1) then - ! -- vz uniform, nonzero + ! vz uniform, nonzero z = zi + vzi * dt else - ! -- vz nonuniform + ! vz nonuniform z = zbot + (vzi * dexp(az * dt) - vzbot) / az end if end if diff --git a/src/Solution/ParticleTracker/TernarySolveTrack.f90 b/src/Solution/ParticleTracker/TernarySolveTrack.f90 index 8a086fccaec..86a418e60ef 100644 --- a/src/Solution/ParticleTracker/TernarySolveTrack.f90 +++ b/src/Solution/ParticleTracker/TernarySolveTrack.f90 @@ -304,14 +304,6 @@ subroutine step_analytical(t, alp, bet) real(DP) :: bet if (icase .eq. 1) then - if (kper > 86) then - print *, "----" - print *, ca1 - print *, ca2 - print *, ca3 - print *, waa - print *, t - end if alp = ca1 + ca2 * dexp(waa * t) + ca3 * dexp(wbb * t) bet = cb1 + cb2 * dexp(wbb * t) else if (icase .eq. -1) then