Skip to content

Commit

Permalink
cleanup, todos
Browse files Browse the repository at this point in the history
  • Loading branch information
@wpbonelli
wpbonelli committed Dec 13, 2023
1 parent 98eba93 commit eac7285
Show file tree
Hide file tree
Showing 5 changed files with 124 additions and 122 deletions.
2 changes: 2 additions & 0 deletions src/Model/ParticleTracking/prt1prp1.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,6 +469,8 @@ subroutine prp_ad(this)
call store_error_unit(this%inunit, terminate=.true.)
end if

! -- todo: make sure release point isn't above or below grid

! -- set particle boundname
if (size(this%boundname) /= 0) &
bndName = this%boundname(nps)
Expand Down
2 changes: 1 addition & 1 deletion src/Solution/ParticleTracker/MethodCellTernary.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ subroutine load_subcell(this, particle, levelNext, subcellTri)
isc = particle%iTrackingDomain(3)
npolyverts = this%cellPoly%defn%npolyverts

! -- Find subcell if not known ! kluge note: from "find_init_triangle"
! -- 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
Expand Down
21 changes: 12 additions & 9 deletions src/Solution/ParticleTracker/MethodSubcellTernary.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,10 +125,12 @@ subroutine track_sub(this, subcell, particle, tmax)

! -- Translate and rotate coordinates to "canonical" configuration
call canonical(x0, y0, x1, y1, x2, y2, &
v0x, v0y, v1x, v1y, v2x, v2y, xi, yi, &
rxx, rxy, ryx, ryy, sxx, sxy, syy, &
lbary, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
v0x, v0y, v1x, v1y, v2x, v2y, &
xi, yi, &
rxx, rxy, ryx, ryy, &
sxx, sxy, syy, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
lbary)

! -- Do calculations related to analytical z solution, which can be done
! -- after traverse_triangle call if results not needed for adaptive time
Expand All @@ -150,11 +152,12 @@ subroutine track_sub(this, subcell, particle, tmax)
! once the final trajectory time is known
call traverse_triangle(ntmax, nsave, diff, rdiff, &
isolv, tol, step, &
dtexitxy, alpexit, betexit, itrifaceenter, &
itrifaceexit, rxx, rxy, ryx, ryy, &
sxx, sxy, syy, lbary, alp0, &
bet0, alp1, bet1, alp2, bet2, &
alpi, beti, vziodz, az)
dtexitxy, alpexit, betexit, &
itrifaceenter, itrifaceexit, &
rxx, rxy, ryx, ryy, &
sxx, sxy, syy, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
vziodz, az, lbary)

! -- Check for no exit face
if ((itopbotexit .eq. 0) .and. (itrifaceexit .eq. 0)) then
Expand Down
209 changes: 101 additions & 108 deletions src/Solution/ParticleTracker/TernarySolveTrack.f90
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
module TernarySolveTrack

use KindModule, only: I4B, DP, LGP
use GeomUtilModule, only: skew
use ErrorUtilModule, only: pstop
private
Expand Down Expand Up @@ -27,17 +28,26 @@ module TernarySolveTrack
!> @brief Traverse triangular cell
subroutine traverse_triangle(ntmax, nsave, diff, rdiff, &
isolv, tol, step, texit, &
alpexit, betexit, itrifaceenter, itrifaceexit, &
rxx, rxy, ryx, ryy, sxx, sxy, syy, &
lbary, alp0, bet0, alp1, bet1, &
alp2, bet2, alpi, beti, vziodz, az)
alpexit, betexit, &
itrifaceenter, itrifaceexit, &
rxx, rxy, ryx, ryy, &
sxx, sxy, syy, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
vziodz, az, &
bary)
use Ternary, only: waa, wab, wba, wbb
! real(DP), intent(inout) :: pexit(2) !< ??
! real(DP), intent(inout) :: itrifaceenter, itrifaceexit !< ??
! real(DP), intent(inout) :: r(2, 2) !< rotation matrix
! real(DP), intent(inout) :: s(2, 2) !< skew matrix
! real(DP), intent(inout) :: t(4, 2) !< triangle points
! logical(LGP), intent(in), optional :: bary !< barycentric coordinates
implicit double precision(a - h, o - z)
logical lbary
logical :: bary
intrinsic cpu_time

! -- Compute elements of matrix W
call get_w(alp0, bet0, alp1, bet1, alp2, bet2, lbary, waa, wab, wba, wbb)
call get_w(alp1, bet1, alp2, bet2, waa, wab, wba, wbb, bary)

! -- Determine alpha and beta analytically as functions of time
call solve_coefs(alpi, beti)
Expand Down Expand Up @@ -76,55 +86,7 @@ subroutine traverse_triangle(ntmax, nsave, diff, rdiff, &

end subroutine

! !> @brief Find initial cell and layer
! !!
! !! kluge: move to "ternary_setup"???
! !!
! !<
! subroutine find_init_cell(numvermx,nvert,nlay,ncpl,xi,yi,zi,ilayeri,icelli)
! !
! use Ternary, only: z_bot
! implicit double precision(a - h, o - z)
! logical point_in_polygon
! !
! ! -- Find initial cell
! icelli = 0
! do icell = 1, ncpl
! if (point_in_polygon(numvermx, nvert, ncpl, xi, yi, icell)) then
! icelli = icell
! exit
! end if
! end do
! if (icelli .eq. 0) then
! write (*, '(A)') "error -- initial cell not found" ! kluge
! write (69, '(A)') "error -- initial cell not found"
! !!pause
! stop
! end if
! !
! ! -- Find initial layer
! ilayeri = 0
! if (zi .gt. z_bot(0, icelli)) then
! write (*, '(A)') "error -- initial z above top of model" ! kluge
! write (69, '(A)') "error -- initial z above top of model"
! !!pause
! stop
! end if
! do ilayer = 1, nlay
! if (zi .ge. z_bot(ilayer, icelli)) then
! ilayeri = ilayer
! exit
! end if
! end do
! if (ilayeri .eq. 0) then
! write (*, '(A)') "error -- initial z below bottom of model" ! kluge
! write (69, '(A)') "error -- initial z below bottom of model"
! !!pause
! stop
! end if
!
! end subroutine

! todo: reinstate
! !> @brief Find initial triangular subcell
! !!
! !! kluge note: not currently used but maybe could be in load_subcell of MethodCellTernary???
Expand Down Expand Up @@ -175,15 +137,31 @@ subroutine traverse_triangle(ntmax, nsave, diff, rdiff, &
!
! end subroutine

!> @brief Translate and rotate coordinates to "canonical" configuration
!> @brief Set coordinates to "canonical" configuration
subroutine canonical(x0, y0, x1, y1, x2, y2, &
v0x, v0y, v1x, v1y, v2x, v2y, &
xi, yi, rxx, rxy, ryx, ryy, sxx, sxy, syy, lbary, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
xi, yi, &
rxx, rxy, ryx, ryy, &
sxx, sxy, syy, &
alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
bary)
use Ternary, only: cv0, cv1, cv2
! real(DP), intent(inout) :: p(3, 2) !< ??
! real(DP), intent(inout) :: v(3, 2) !< ??
! real(DP), intent(inout) :: i(2) !< ??
! real(DP), intent(inout) :: r(2, 2) !< rotation matrix
! real(DP), intent(inout) :: s(2, 2) !< skew matrix
! real(DP), intent(inout) :: t(3, 2) !< triangle points
! logical(LGP), intent(in), optional :: bary !< barycentric coordinates
implicit double precision(a - h, o - z)
logical lbary
double precision :: rot(2, 2), res(2)
logical :: bary
real(DP) :: rot(2, 2), res(2)

! if (present(bary)) then
! lbary = bary
! else
! lbary = .true.
! end if

! -- Translate and rotate coordinates to "canonical" configuration
x1diff = x1 - x0
Expand All @@ -203,6 +181,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, &
alp1 = baselen
bet1 = 0d0

! todo refactor alp/bet vars
rot = reshape((/rxx, ryx, rxy, ryy/), shape(rot))
res = matmul(rot, (/x2diff, y2diff/))
alp2 = res(1)
Expand All @@ -218,7 +197,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, &
alpi = res(1)
beti = res(2)

if (lbary) then
if (bary) then
sxx = 1d0 / alp1
syy = 1d0 / bet2
sxy = -alp2 * sxx * syy
Expand All @@ -228,6 +207,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, &
call skew(sxx, sxy, syy, cv0)
call skew(sxx, sxy, syy, cv1)
call skew(sxx, sxy, syy, cv2)
! todo turn alpi/beti into array
res = (/alpi, beti/)
call skew(sxx, sxy, syy, res)
alpi = res(1)
Expand All @@ -237,17 +217,34 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, &
end subroutine

!> @brief Compute elements of W matrix
subroutine get_w(alp0, bet0, alp1, bet1, alp2, bet2, lbary, waa, wab, wba, wbb)
subroutine get_w( &
alp1, bet1, alp2, bet2, &
waa, wab, wba, wbb, &
bary)
use Ternary, only: cv0, cv1, cv2
implicit double precision(a - h, o - z)
logical lbary
! -- dummy
real(DP), intent(in) :: alp1, bet1, alp2, bet2 !< triangle face points
real(DP), intent(inout) :: waa, wab, wba, wbb !< w matrix
! real(DP), intent(inout) :: t1(2), t2(2) !< triangle face points
! real(DP), intent(inout) :: w(2, 2) !< w matrix
logical(LGP), intent(in), optional :: bary !< barycentric coordinates
! -- local
logical(LGP) :: lbary
real(DP) :: v1alpdiff, v2alpdiff, v2betdiff
real(DP) :: ooalp1, oobet2, vterm

if (present(bary)) then
lbary = bary
else
lbary = .true.
end if

! -- Note: wab is the "alpha,beta" entry in matrix W
! and the alpha component of the w^(beta) vector
v1alpdiff = cv1(1) - cv0(1)
v2alpdiff = cv2(1) - cv0(1)
v2betdiff = cv2(2) - cv0(2)
if (lbary) then
if (bary) then
waa = v1alpdiff
wab = v2alpdiff
wba = 0d0
Expand Down Expand Up @@ -839,7 +836,7 @@ subroutine soln_euler(itriface, alpi, beti, step, vziodz, &
!! Simulation of Classical and Quantum Systems," 2nd ed., Springer, New York.
!!
!<
double precision function zeroch(x0, x1, f, epsa)
function zeroch(x0, x1, f, epsa)
implicit double precision(a - h, o - z)

epsm = epsilon(x0)
Expand Down Expand Up @@ -924,14 +921,9 @@ double precision function zeroch(x0, x1, f, epsa)
end function

!> @brief Compute a zero of the function f(x) in the interval (x0, x1)
!!
!! A zero of the function f{x} is computed in the interval (x0, x1)
!! given tolerance epsa using a test method.
!!
!<
double precision function zerotest(x0, x1, f, epsa) ! kluge test
function zerotest(x0, x1, f, epsa)
implicit double precision(a - h, o - z)
logical retainedxa, retainedxb
logical(LGP) :: retainedxa, retainedxb

epsm = epsilon(x0)
f0 = f(x0)
Expand Down Expand Up @@ -998,51 +990,54 @@ double precision function zerotest(x0, x1, f, epsa) ! kluge test

end function

!> @brief Compute a zero of the function f(x) in the interval (x0, x1)
double precision function zeroin(ax, bx, f, tol)
double precision ax, bx, f, tol
! a zero of the function f(x) is computed in the interval ax,bx .
! input..
!
! ax left endpoint of initial interval
! bx right endpoint of initial interval
! f function subprogram which evaluates f(x) for any x in
! the interval ax,bx
! tol desired length of the interval of uncertainty of the
! final result (.ge.0.)
!
! output..
!
! zeroin abscissa approximating a zero of f in the interval ax,bx
!
! it is assumed that f(ax) and f(bx) have opposite signs
! this is checked, and an error message is printed if this is not
! satisfied. zeroin returns a zero x in the given interval
! ax,bx to within a tolerance 4*macheps*abs(x)+tol, where macheps is
! the relative machine precision defined as the smallest representable
! number such that 1.+macheps .gt. 1.
! this function subprogram is a slightly modified translation of
! the algol 60 procedure zero given in richard brent, algorithms for
! minimization without derivatives, prentice-hall, inc. (1973).
double precision a, b, c, d, e, eps, fa, fb, fc, tol1, xm, p, q, r, s
!> @brief Compute a zero of the function f(x) in the interval (x0, x1).
!!
!! A zero of the function f(x) is computed in the interval ax,bx .
!!
!! input..
!!
!! ax left endpoint of initial interval
!! bx right endpoint of initial interval
!! f function subprogram which evaluates f(x) for any x in
!! the interval ax,bx
!! tol desired length of the interval of uncertainty of the
!! final result (.ge.0.)
!!
!! output..
!!
!! zeroin abscissa approximating a zero of f in the interval ax,bx
!!
!! it is assumed that f(ax) and f(bx) have opposite signs
!! this is checked, and an error message is printed if this is not
!! satisfied. zeroin returns a zero x in the given interval
!! ax,bx to within a tolerance 4*macheps*abs(x)+tol, where macheps is
!! the relative machine precision defined as the smallest representable
!! number such that 1.+macheps .gt. 1.
!! this function subprogram is a slightly modified translation of
!! the algol 60 procedure zero given in richard brent, algorithms for
!! minimization without derivatives, prentice-hall, inc. (1973).
!<
function zeroin(ax, bx, f, tol)
implicit double precision(a - h, o - z)

eps = epsilon(ax)
tol1 = eps + 1.0d0

a = ax
b = bx
fa = f(a)
fb = f(b)

! check that f(ax) and f(bx) have different signs
if (fa .eq. 0.0d0 .or. fb .eq. 0.0d0) go to 20
if (fa * (fb / dabs(fb)) .le. 0.0d0) go to 20
write (6, 2500)
2500 format(1x, 'f(ax) and f(bx) do not have different signs,', &
' zeroin is aborting')
return
if (.not. ((fa .eq. 0.0d0 .or. fb .eq. 0.0d0) .or. &
(fa * (fb / dabs(fb)) .le. 0.0d0))) &
call pstop(1, 'f(ax) and f(bx) do not have different signs,')

20 c = a
fc = fa
d = b - a
e = d

30 if (dabs(fc) .ge. dabs(fb)) go to 40
a = b
b = c
Expand All @@ -1055,7 +1050,6 @@ double precision function zeroin(ax, bx, f, tol)
if ((dabs(xm) .le. tol1) .or. (fb .eq. 0.0d0)) go to 150

! see if a bisection is forced

if ((dabs(e) .ge. tol1) .and. (dabs(fa) .gt. dabs(fb))) go to 50
d = xm
e = d
Expand All @@ -1068,7 +1062,6 @@ double precision function zeroin(ax, bx, f, tol)
go to 70

! inverse quadratic interpolation

60 q = fa / fc
r = fb / fc
p = s * (2.0d0 * xm * q * (q - r) - (b - a) * (r - 1.0d0))
Expand Down
Loading

0 comments on commit eac7285

Please sign in to comment.