fix(prt): convert to local coords in ternary cell method (#2122)

This reverts the switch to calculating the centroid via arithmetic mean for triangular cells in #2119 — the rest of the vertex velocity calculation assumes the polygonal definition of the centroid so that change was an error. A more robust way to avoid intermediate loss of precision is converting to local cell-local coordinates before applying the ternary cell method.

Author: wpbonelli

src/Solution/ParticleTracker/MethodCellTernary.f90

@@ -8,7 +8,7 @@ module MethodCellTernaryModule
   use CellDefnModule
   use SubcellTriModule, only: SubcellTriType, create_subcell_tri
   use ParticleModule
-  use GeomUtilModule, only: area
+  use GeomUtilModule, only: area, transform
   use ConstantsModule, only: DZERO, DONE, DTWO, DSIX
   use GeomUtilModule, only: point_in_polygon
   implicit none
@@ -158,6 +158,8 @@ subroutine apply_mct(this, particle, tmax)
     real(DP), intent(in) :: tmax
     ! local
     integer(I4B) :: i
+    real(DP) :: x, y, z, xO, yO
+    real(DP), allocatable :: xs(:), ys(:)
 
     ! (Re)allocate type-bound arrays
     select type (cell => this%cell)
@@ -194,10 +196,24 @@ subroutine apply_mct(this, particle, tmax)
       allocate (this%xvertnext(this%nverts))
       allocate (this%yvertnext(this%nverts))
 
+      ! New origin is the corner of the cell's
+      ! bounding box closest to the old origin
+      allocate (xs(this%nverts))
+      allocate (ys(this%nverts))
+      xs = cell%defn%polyvert(1, :)
+      ys = cell%defn%polyvert(2, :)
+      xO = xs(minloc(abs(xs), dim=1))
+      yO = ys(minloc(abs(ys), dim=1))
+      deallocate (xs)
+      deallocate (ys)
+
       ! Cell vertices
       do i = 1, this%nverts
-        this%xvert(i) = cell%defn%polyvert(1, i)
-        this%yvert(i) = cell%defn%polyvert(2, i)
+        x = cell%defn%polyvert(1, i)
+        y = cell%defn%polyvert(2, i)
+        call transform(x, y, DZERO, x, y, z, xO, yO)
+        this%xvert(i) = x
+        this%yvert(i) = y
       end do
 
       ! Top, bottom, and thickness
@@ -212,13 +228,21 @@ subroutine apply_mct(this, particle, tmax)
       this%iprev = cshift(this%iprev, -1)
       this%xvertnext = cshift(this%xvert, 1)
       this%yvertnext = cshift(this%yvert, 1)
-    end select
 
-    ! Calculate vertex velocities.
-    call this%vertvelo()
+      ! Calculate vertex velocities.
+      call this%vertvelo()
+
+      ! Transform particle coordinates
+      call particle%transform(xO, yO)
+
+      ! Track the particle across the cell.
+      call this%track(particle, 2, tmax)
 
-    ! Track the particle across the cell.
-    call this%track(particle, 2, tmax)
+      ! Transform particle coordinates back
+      call particle%transform(xO, yO, invert=.true.)
+      call particle%reset_transform()
+
+    end select
 
   end subroutine apply_mct
 
@@ -403,16 +427,11 @@ subroutine vertvelo(this)
       sixa = areacell * DSIX
       wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert)
       nvert = size(this%xvert)
-      if (nvert == 3) then
-        this%xctr = sum(this%xvert) / 3.0_DP
-        this%yctr = sum(this%yvert) / 3.0_DP
-      else
-        this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
-        this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
-      end if
+      this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
+      this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
 
-      ! TODO: do we need to check if center is in polygon?
-      ! only if using the centroid method which misbehaves
+      ! TODO: can we use some of the terms of the centroid calculation
+      ! to do a cheap point in polygon check?
       !
       ! allocate(poly(2, nvert))
       ! poly(1,:) = this%xvert