From eac7285b2c3bae3628a7f406161fa87cc6db1a0d Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Tue, 12 Dec 2023 19:56:16 -0500 Subject: [PATCH] cleanup, todos --- src/Model/ParticleTracking/prt1prp1.f90 | 2 + .../ParticleTracker/MethodCellTernary.f90 | 2 +- .../ParticleTracker/MethodSubcellTernary.f90 | 21 +- .../ParticleTracker/TernarySolveTrack.f90 | 209 +++++++++--------- src/Utilities/GeomUtil.f90 | 12 +- 5 files changed, 124 insertions(+), 122 deletions(-) diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index 44938f39a42..89e869ed327 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -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) diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index a32f86da94a..386e3d3c658 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -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 diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index a965cd96543..94ea3f39bc9 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -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 @@ -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 diff --git a/src/Solution/ParticleTracker/TernarySolveTrack.f90 b/src/Solution/ParticleTracker/TernarySolveTrack.f90 index cafac0c340a..4cb06c59614 100644 --- a/src/Solution/ParticleTracker/TernarySolveTrack.f90 +++ b/src/Solution/ParticleTracker/TernarySolveTrack.f90 @@ -1,5 +1,6 @@ module TernarySolveTrack + use KindModule, only: I4B, DP, LGP use GeomUtilModule, only: skew use ErrorUtilModule, only: pstop private @@ -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) @@ -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??? @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 @@ -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) @@ -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) @@ -998,33 +990,36 @@ 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 @@ -1032,17 +1027,17 @@ double precision function zeroin(ax, bx, f, tol) 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 @@ -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 @@ -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)) diff --git a/src/Utilities/GeomUtil.f90 b/src/Utilities/GeomUtil.f90 index e2b98604836..40c2e19b347 100644 --- a/src/Utilities/GeomUtil.f90 +++ b/src/Utilities/GeomUtil.f90 @@ -179,9 +179,11 @@ subroutine skew(sxx, sxy, syy, vec, invert) end subroutine skew !> @brief Apply 3D translation and optional 2D rotation to coordinates. - subroutine transform_coords(xin, yin, zin, xout, yout, zout, & + subroutine transform_coords(xin, yin, zin, & + xout, yout, zout, & xorigin, yorigin, zorigin, & - sinrot, cosrot, invert) + sinrot, cosrot, & + invert) ! -- dummy real(DP) :: xin, yin, zin !< input coordinates real(DP) :: xout, yout, zout !< output coordinates @@ -237,9 +239,11 @@ subroutine transform_coords(xin, yin, zin, xout, yout, zout, & end subroutine transform_coords !> @brief Apply a 3D translation and 2D rotation to an existing transformation. - subroutine modify_transf(xorigin, yorigin, zorigin, sinrot, cosrot, & + subroutine modify_transf(xorigin, yorigin, zorigin, & + sinrot, cosrot, & xorigin_new, yorigin_new, zorigin_new, & - sinrot_new, cosrot_new, invert) + sinrot_new, cosrot_new, & + invert) ! -- dummy real(DP) :: xorigin, yorigin, zorigin !< origin coordinates (original) real(DP) :: sinrot, cosrot !< sine and cosine of rotation (original)