coord_transform/m_coord_transform_frenet_serret.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Frenet-Serret coordinate transformation | ||
| 2 | !> | ||
| 3 | !> This module implements a coordinate transformation based on a Frenet-Serret frame | ||
| 4 | !> along a space curve γ(zp). The transformation is defined as: | ||
| 5 | !> | ||
| 6 | !> x(xp, yp, zp) = γ(zp) - R(xp, yp) * S(yp) * N(zp) + R(xp, yp) * C(yp) * M(zp) | ||
| 7 | !> | ||
| 8 | !> where: | ||
| 9 | !> - γ(zp) is the base curve parameterized by zp | ||
| 10 | !> - T(zp), N(zp), M(zp) are the tangent, normal, and binormal vectors (Frenet-Serret frame) | ||
| 11 | !> - R(xp, yp), C(yp), and S(yp) are functions defining the cross-sectional shape | ||
| 12 | !> | ||
| 13 | !> The Frenet-Serret frame is computed from the curve γ using: | ||
| 14 | !> - T = γ' / ||γ'|| (unit tangent) | ||
| 15 | !> - N = T' / ||T'|| (unit normal) | ||
| 16 | !> - M = T × N (unit binormal) | ||
| 17 | !> - κ = ||γ' × γ''|| / ||γ'||³ (curvature) | ||
| 18 | !> - τ = (γ' × γ'') · γ''' / ||γ' × γ''||² (torsion) | ||
| 19 | !> | ||
| 20 | !> Parameters: | ||
| 21 | !> - R0_pert :: real(wp) | ||
| 22 | !> - R1 :: real(wp) | ||
| 23 | !> - R1_pert :: real(wp) | ||
| 24 | !> - n0 :: integer | ||
| 25 | !> - n1 :: integer | ||
| 26 | !> | ||
| 27 | !> @note This module was automatically generated by frenet_serret_generate.jl | ||
| 28 | module m_coord_transform_frenet_serret | ||
| 29 | use m_common, only: wp, pi | ||
| 30 | use m_coord_transform_abstract, only: CoordTransformAbstract | ||
| 31 | |||
| 32 | implicit none | ||
| 33 | |||
| 34 | private | ||
| 35 | public :: FrenetSerretTransform | ||
| 36 | public :: frenet_serret_transform, frenet_serret_inverse_transform, frenet_serret_jacobian | ||
| 37 | public :: frenet_serret_jacobian_matrix, frenet_serret_jacobian_matrix_inv | ||
| 38 | public :: frenet_serret_G_matrix, frenet_serret_G_matrix_inv | ||
| 39 | |||
| 40 | !> @brief Derived type for Frenet-Serret coordinate transformation | ||
| 41 | type, extends(CoordTransformAbstract) :: FrenetSerretTransform | ||
| 42 | real(wp) :: R0_pert !< Parameter from the coordinate transformation | ||
| 43 | real(wp) :: R1 !< Parameter from the coordinate transformation | ||
| 44 | real(wp) :: R1_pert !< Parameter from the coordinate transformation | ||
| 45 | integer :: n0 !< Parameter from the coordinate transformation | ||
| 46 | integer :: n1 !< Parameter from the coordinate transformation | ||
| 47 | contains | ||
| 48 | ! Additional helpers for efficiency | ||
| 49 | procedure :: curvature => frenet_serret_curvature_tbp | ||
| 50 | procedure :: torsion => frenet_serret_torsion_tbp | ||
| 51 | procedure :: tangent => frenet_serret_tangent_tbp | ||
| 52 | procedure :: normal => frenet_serret_normal_tbp | ||
| 53 | procedure :: binormal => frenet_serret_binormal_tbp | ||
| 54 | procedure :: derivative_magnitude => frenet_serret_derivative_magnitude_tbp | ||
| 55 | end type FrenetSerretTransform | ||
| 56 | |||
| 57 | interface FrenetSerretTransform | ||
| 58 | module procedure init_FrenetSerretTransform | ||
| 59 | end interface FrenetSerretTransform | ||
| 60 | |||
| 61 | contains | ||
| 62 | |||
| 63 | 7 | function init_FrenetSerretTransform(R0_pert, R1, R1_pert, n0, n1) result(this) | |
| 64 | implicit none | ||
| 65 | real(wp), intent(in) :: R0_pert | ||
| 66 | real(wp), intent(in) :: R1 | ||
| 67 | real(wp), intent(in) :: R1_pert | ||
| 68 | integer, intent(in) :: n0 | ||
| 69 | integer, intent(in) :: n1 | ||
| 70 | type(FrenetSerretTransform) :: this | ||
| 71 | |||
| 72 | this%is_orthogonal = .false. | ||
| 73 | this%has_polar_xy_singularity = .true. | ||
| 74 | 7 | this%R0_pert = R0_pert | |
| 75 | 7 | this%R1 = R1 | |
| 76 | 7 | this%R1_pert = R1_pert | |
| 77 | 7 | this%n0 = n0 | |
| 78 | 7 | this%n1 = n1 | |
| 79 | 7 | end function init_FrenetSerretTransform | |
| 80 | |||
| 81 | 2100 | pure real(wp) function frenet_serret_transform(this, i, xp, yp, zp) result(ans) | |
| 82 | !$acc routine seq | ||
| 83 | implicit none | ||
| 84 | type(FrenetSerretTransform), intent(in) :: this | ||
| 85 | integer, intent(in) :: i | ||
| 86 | real(wp), intent(in) :: xp, yp, zp | ||
| 87 | real(wp) :: gamma_val, N_val, M_val, Rval, Cval, Sval | ||
| 88 | |||
| 89 | ! Use helper TBPs for gamma, R, C, S, normal, and binormal | ||
| 90 | 2100 | gamma_val = frenet_serret_curve(this, i, zp) | |
| 91 | 2100 | Rval = frenet_serret_R(this, xp, yp, zp) | |
| 92 | Cval = frenet_serret_C(this, xp, yp, zp) | ||
| 93 | Sval = frenet_serret_S(this, xp, yp, zp) | ||
| 94 | 2100 | N_val = frenet_serret_normal(this, i, zp) | |
| 95 | 2100 | M_val = frenet_serret_binormal(this, i, zp) | |
| 96 | |||
| 97 | 2100 | ans = gamma_val - Rval * Sval * N_val + Rval * Cval * M_val | |
| 98 | 2100 | end function frenet_serret_transform | |
| 99 | |||
| 100 | 7200 | pure real(wp) function frenet_serret_jacobian_matrix(this, i, j, xp, yp, zp) result(ans) | |
| 101 | !$acc routine seq | ||
| 102 | implicit none | ||
| 103 | type(FrenetSerretTransform), intent(in) :: this | ||
| 104 | integer, intent(in) :: i, j | ||
| 105 | real(wp), intent(in) :: xp, yp, zp | ||
| 106 | real(wp) :: Rval, Cval, Sval, Rx_val, Ry_val, Cy_val, Sy_val | ||
| 107 | real(wp) :: T_i, N_i, M_i, gp_i, kappa_val, tau_val, gpnorm_val | ||
| 108 | |||
| 109 | ! Compute R, C, S and their derivatives using helper TBPs | ||
| 110 | 7200 | Rval = frenet_serret_R(this, xp, yp, zp) | |
| 111 | Cval = frenet_serret_C(this, xp, yp, zp) | ||
| 112 | Sval = frenet_serret_S(this, xp, yp, zp) | ||
| 113 | 7200 | Rx_val = frenet_serret_Rx(this, xp, yp, zp) | |
| 114 | 7200 | Ry_val = frenet_serret_Ry(this, xp, yp, zp) | |
| 115 | Cy_val = frenet_serret_Cy(this, xp, yp, zp) | ||
| 116 | Sy_val = frenet_serret_Sy(this, xp, yp, zp) | ||
| 117 | |||
| 118 | ! Derivatives: ∂(R*S)/∂x = Rx*S; ∂(R*S)/∂y = Ry*S + R*Sy | ||
| 119 | ! ∂(R*C)/∂x = Rx*C; ∂(R*C)/∂y = Ry*C + R*Cy | ||
| 120 | |||
| 121 | ! Get Frenet-Serret frame components and geometric properties | ||
| 122 | 7200 | T_i = frenet_serret_tangent(this, i, zp) | |
| 123 | 7200 | N_i = frenet_serret_normal(this, i, zp) | |
| 124 | 7200 | M_i = frenet_serret_binormal(this, i, zp) | |
| 125 | 7200 | kappa_val = frenet_serret_curvature(this, zp) | |
| 126 | 7200 | tau_val = frenet_serret_torsion(this, zp) | |
| 127 | 7200 | gpnorm_val = frenet_serret_derivative_magnitude(this, zp) | |
| 128 | |||
| 129 | 9600 | select case (j) | |
| 130 | case (1) ! ∂r/∂x = -∂R/∂x S N + ∂R/∂x C M | ||
| 131 | 2400 | ans = -(Rx_val * Sval) * N_i + (Rx_val * Cval) * M_i | |
| 132 | case (2) ! ∂r/∂y = -∂(R*S)/∂y N + ∂(R*C)/∂y M | ||
| 133 | 2400 | ans = -(Ry_val * Sval + Rval * Sy_val) * N_i + (Ry_val * Cval + Rval * Cy_val) * M_i | |
| 134 | case (3) ! ∂r/∂z = γ' - (R*S) dN/dz + (R*C) dM/dz | ||
| 135 | ! Using Frenet-Serret: dN/dz = (-κT + τM)||γ'||, dM/dz = -τN||γ'|| | ||
| 136 | 2400 | gp_i = frenet_serret_curve_derivative(this, i, zp) | |
| 137 | 2400 | ans = gp_i - (Rval * Sval) * (-kappa_val * T_i + tau_val * M_i) * gpnorm_val + (Rval * Cval) * (-tau_val * N_i) * gpnorm_val | |
| 138 | end select | ||
| 139 | 7200 | end function frenet_serret_jacobian_matrix | |
| 140 | |||
| 141 | 22132424 | pure real(wp) function frenet_serret_jacobian(this, xp, yp, zp) result(ans) | |
| 142 | !$acc routine seq | ||
| 143 | implicit none | ||
| 144 | type(FrenetSerretTransform), intent(in) :: this | ||
| 145 | real(wp), intent(in) :: xp, yp, zp | ||
| 146 | real(wp) :: Rval, Sval, kappa_val, gpnorm_val, det2x2 | ||
| 147 | |||
| 148 | ! Use the relation: detJ = detJ2D * ||γ'|| * (1 + κ*R*S) | ||
| 149 | 22132424 | det2x2 = frenet_serret_det2D(this, xp, yp, zp) | |
| 150 | 22132424 | gpnorm_val = frenet_serret_derivative_magnitude(this, zp) | |
| 151 | 22132424 | kappa_val = frenet_serret_curvature(this, zp) | |
| 152 | 22132424 | Rval = frenet_serret_R(this, xp, yp, zp) | |
| 153 | Sval = frenet_serret_S(this, xp, yp, zp) | ||
| 154 | |||
| 155 | 22132424 | ans = det2x2 * gpnorm_val * (1.0_wp + kappa_val * Rval * Sval) | |
| 156 | 22132424 | end function frenet_serret_jacobian | |
| 157 | |||
| 158 | 900 | pure real(wp) function frenet_serret_jacobian_matrix_inv(this, i, j, xp, yp, zp) result(ans) | |
| 159 | !$acc routine seq | ||
| 160 | implicit none | ||
| 161 | type(FrenetSerretTransform), intent(in) :: this | ||
| 162 | integer, intent(in) :: i, j | ||
| 163 | real(wp), intent(in) :: xp, yp, zp | ||
| 164 | integer :: k | ||
| 165 | |||
| 166 | ! Use the identity: J^(-1) = G^(-1) * J^T | ||
| 167 | ! So: Jinv(i,j) = sum_k Ginv(i,k) * J(j,k) | ||
| 168 | |||
| 169 | ! For non-orthogonal coordinates: | ||
| 170 | ! G has structure with G(1,3) = 0 but G(2,3) ≠ 0 | ||
| 171 | ! Use full 3x3 sum: Jinv(i,j) = sum_k Ginv(i,k) * J(j,k) | ||
| 172 | |||
| 173 | ans = 0._wp | ||
| 174 |
2/2✓ Branch 3 → 4 taken 2700 times.
✓ Branch 3 → 5 taken 900 times.
|
3600 | do k = 1, 3 |
| 175 | 3600 | ans = ans + frenet_serret_G_matrix_inv(this, i, k, xp, yp, zp) * frenet_serret_jacobian_matrix(this, j, k, xp, yp, zp) | |
| 176 | end do | ||
| 177 | 900 | end function frenet_serret_jacobian_matrix_inv | |
| 178 | |||
| 179 | 146692440 | pure real(wp) function frenet_serret_G_matrix(this, i, j, xp, yp, zp) result(ans) | |
| 180 | !$acc routine seq | ||
| 181 | implicit none | ||
| 182 | type(FrenetSerretTransform), intent(in) :: this | ||
| 183 | integer, intent(in) :: i, j | ||
| 184 | real(wp), intent(in) :: xp, yp, zp | ||
| 185 | real(wp) :: Rval, Cval, Sval, Rx_val, Ry_val, Cy_val, Sy_val | ||
| 186 | real(wp) :: kappa_val, gpnorm_val, tau_val | ||
| 187 | |||
| 188 | ! Compute R, C, S and their derivatives using helper TBPs | ||
| 189 | 146692440 | Rval = frenet_serret_R(this, xp, yp, zp) | |
| 190 | Cval = frenet_serret_C(this, xp, yp, zp) | ||
| 191 | Sval = frenet_serret_S(this, xp, yp, zp) | ||
| 192 | 146692440 | Rx_val = frenet_serret_Rx(this, xp, yp, zp) | |
| 193 | 146692440 | Ry_val = frenet_serret_Ry(this, xp, yp, zp) | |
| 194 | Cy_val = frenet_serret_Cy(this, xp, yp, zp) | ||
| 195 | Sy_val = frenet_serret_Sy(this, xp, yp, zp) | ||
| 196 | |||
| 197 | ! Get geometric properties | ||
| 198 | 146692440 | kappa_val = frenet_serret_curvature(this, zp) | |
| 199 | 146692440 | tau_val = frenet_serret_torsion(this, zp) | |
| 200 | 146692440 | gpnorm_val = frenet_serret_derivative_magnitude(this, zp) | |
| 201 | |||
| 202 | 146692440 | if (i == j) then | |
| 203 | 29338328 | select case (i) | |
| 204 | case (1) | ||
| 205 | 29338328 | ans = Rx_val**2 * (Sval**2 + Cval**2) | |
| 206 | case (2) | ||
| 207 | 29338328 | ans = (Ry_val * Sval + Rval * Sy_val)**2 + (Ry_val * Cval + Rval * Cy_val)**2 | |
| 208 | case (3) | ||
| 209 |
3/4✓ Branch 3 → 4 taken 29338328 times.
✓ Branch 3 → 5 taken 29338328 times.
✓ Branch 3 → 6 taken 29338328 times.
✗ Branch 3 → 21 not taken.
|
88014984 | ans = gpnorm_val**2 * ((1.0_wp + kappa_val * (Rval * Sval))**2 + (Rval * tau_val)**2 * (Cval**2 + Sval**2)) |
| 210 | end select | ||
| 211 | else | ||
| 212 |
8/8✓ Branch 7 → 8 taken 29338528 times.
✓ Branch 7 → 9 taken 29338928 times.
✓ Branch 8 → 9 taken 200 times.
✓ Branch 8 → 11 taken 29338328 times.
✓ Branch 9 → 10 taken 29338528 times.
✓ Branch 9 → 12 taken 600 times.
✓ Branch 10 → 11 taken 200 times.
✓ Branch 10 → 12 taken 29338328 times.
|
58677456 | if ((i == 1 .and. j == 2) .or. (i == 2 .and. j == 1)) then |
| 213 | 29338528 | ans = Rx_val * Ry_val * (Sval**2 + Cval**2) + Rval * Rx_val * (Cval * Cy_val + Sval * Sy_val) | |
| 214 |
7/8✓ Branch 12 → 13 taken 200 times.
✓ Branch 12 → 14 taken 29338728 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 21 taken 200 times.
✓ Branch 14 → 15 taken 400 times.
✓ Branch 14 → 16 taken 29338328 times.
✓ Branch 15 → 16 taken 200 times.
✓ Branch 15 → 21 taken 200 times.
|
29338928 | else if ((i == 1 .and. j == 3) .or. (i == 3 .and. j == 1)) then |
| 215 | ! G(1,3) = 0 by construction of r = γ - R S N + R C M | ||
| 216 | ans = 0._wp | ||
| 217 |
5/8✓ Branch 16 → 17 taken 29338328 times.
✓ Branch 16 → 18 taken 200 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 20 taken 29338328 times.
✓ Branch 18 → 19 taken 200 times.
✗ Branch 18 → 21 not taken.
✓ Branch 19 → 20 taken 200 times.
✗ Branch 19 → 21 not taken.
|
29338528 | else if ((i == 2 .and. j == 3) .or. (i == 3 .and. j == 2)) then |
| 218 | 29338528 | ans = Rval**2 * gpnorm_val * tau_val * (Cval * Sy_val - Sval * Cy_val) | |
| 219 | end if | ||
| 220 | end if | ||
| 221 | 146692440 | end function frenet_serret_G_matrix | |
| 222 | |||
| 223 | 29338128 | pure real(wp) function frenet_serret_G_matrix_inv(this, i, j, xp, yp, zp) result(ans) | |
| 224 | !$acc routine seq | ||
| 225 | implicit none | ||
| 226 | type(FrenetSerretTransform), intent(in) :: this | ||
| 227 | integer, intent(in) :: i, j | ||
| 228 | real(wp), intent(in) :: xp, yp, zp | ||
| 229 | real(wp) :: G11, G12, G13, G22, G23, G33, det | ||
| 230 | |||
| 231 | ! G has structure: | ||
| 232 | ! G = [ G11 G12 0 ] | ||
| 233 | ! [ G12 G22 G23 ] | ||
| 234 | ! [ 0 G23 G33 ] | ||
| 235 | ! | ||
| 236 | ! The inverse of this matrix requires computing the full determinant | ||
| 237 | ! and cofactor matrix. For efficiency, we compute only the requested component. | ||
| 238 | |||
| 239 | ! Get all non-zero components of G | ||
| 240 | 29338128 | G11 = frenet_serret_G_matrix(this, 1, 1, xp, yp, zp) | |
| 241 | 29338128 | G12 = frenet_serret_G_matrix(this, 1, 2, xp, yp, zp) | |
| 242 | 29338128 | G22 = frenet_serret_G_matrix(this, 2, 2, xp, yp, zp) | |
| 243 | 29338128 | G23 = frenet_serret_G_matrix(this, 2, 3, xp, yp, zp) | |
| 244 | 29338128 | G33 = frenet_serret_G_matrix(this, 3, 3, xp, yp, zp) | |
| 245 | ! Compute determinant: det(G) = G11*(G22*G33 - G23^2) - G12^2*G33 | ||
| 246 | 29338128 | det = G11 * (G22 * G33 - G23**2) - G12**2 * G33 | |
| 247 | |||
| 248 | ! Compute the requested component using cofactors | ||
| 249 |
4/4✓ Branch 2 → 3 taken 9779376 times.
✓ Branch 2 → 5 taken 19558752 times.
✓ Branch 3 → 4 taken 3023248 times.
✓ Branch 3 → 5 taken 6756128 times.
|
29338128 | if (i == 1 .and. j == 1) then |
| 250 | ! Ginv(1,1) = (G22*G33 - G23^2) / det | ||
| 251 | 3023248 | ans = (G22 * G33 - G23**2) / det | |
| 252 |
8/8✓ Branch 5 → 6 taken 6756128 times.
✓ Branch 5 → 7 taken 19558752 times.
✓ Branch 6 → 7 taken 3378064 times.
✓ Branch 6 → 9 taken 3378064 times.
✓ Branch 7 → 8 taken 9779376 times.
✓ Branch 7 → 10 taken 13157440 times.
✓ Branch 8 → 9 taken 3378064 times.
✓ Branch 8 → 10 taken 6401312 times.
|
26314880 | else if ((i == 1 .and. j == 2) .or. (i == 2 .and. j == 1)) then |
| 253 | ! Ginv(1,2) = Ginv(2,1) = -G12*G33 / det | ||
| 254 | 6756128 | ans = -G12 * G33 / det | |
| 255 |
7/8✓ Branch 10 → 11 taken 3378064 times.
✓ Branch 10 → 12 taken 16180688 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 14 taken 3378064 times.
✓ Branch 12 → 13 taken 9779376 times.
✓ Branch 12 → 15 taken 6401312 times.
✓ Branch 13 → 14 taken 3378064 times.
✓ Branch 13 → 15 taken 6401312 times.
|
19558752 | else if ((i == 1 .and. j == 3) .or. (i == 3 .and. j == 1)) then |
| 256 | ! Ginv(1,3) = Ginv(3,1) = G12*G23 / det | ||
| 257 | 6756128 | ans = G12 * G23 / det | |
| 258 |
4/4✓ Branch 15 → 16 taken 6401312 times.
✓ Branch 15 → 18 taken 6401312 times.
✓ Branch 16 → 17 taken 3023248 times.
✓ Branch 16 → 18 taken 3378064 times.
|
12802624 | else if (i == 2 .and. j == 2) then |
| 259 | ! Ginv(2,2) = G11*G33 / det | ||
| 260 | 3023248 | ans = G11 * G33 / det | |
| 261 |
6/8✓ Branch 18 → 19 taken 3378064 times.
✓ Branch 18 → 20 taken 6401312 times.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 22 taken 3378064 times.
✓ Branch 20 → 21 taken 6401312 times.
✗ Branch 20 → 23 not taken.
✓ Branch 21 → 22 taken 3378064 times.
✓ Branch 21 → 23 taken 3023248 times.
|
9779376 | else if ((i == 2 .and. j == 3) .or. (i == 3 .and. j == 2)) then |
| 262 | ! Ginv(2,3) = Ginv(3,2) = -G11*G23 / det | ||
| 263 | 6756128 | ans = -G11 * G23 / det | |
| 264 |
2/4✓ Branch 23 → 24 taken 3023248 times.
✗ Branch 23 → 26 not taken.
✓ Branch 24 → 25 taken 3023248 times.
✗ Branch 24 → 26 not taken.
|
3023248 | else if (i == 3 .and. j == 3) then |
| 265 | ! Ginv(3,3) = (G11*G22 - G12^2) / det | ||
| 266 | 3023248 | ans = (G11 * G22 - G12**2) / det | |
| 267 | end if | ||
| 268 | 29338128 | end function frenet_serret_G_matrix_inv | |
| 269 | |||
| 270 | ! -- curvature and torsion -------------------------------------------- | ||
| 271 | 168832264 | pure real(wp) function frenet_serret_curvature(this, zp) result(ans) | |
| 272 | !$acc routine seq | ||
| 273 | implicit none | ||
| 274 | type(FrenetSerretTransform), intent(in) :: this | ||
| 275 | real(wp), intent(in) :: zp | ||
| 276 | real(wp) :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9 | ||
| 277 | |||
| 278 | 168832264 | tmp1 = this%R0_pert**2.0_wp | |
| 279 | 168832264 | tmp2 = tmp1 * this%R0_pert | |
| 280 | 168832264 | tmp3 = tmp2 * this%R0_pert | |
| 281 | 168832264 | tmp4 = this%n0**2.0_wp | |
| 282 | 168832264 | tmp5 = tmp4 * tmp4 | |
| 283 | 168832264 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 284 | 168832264 | tmp7 = tmp6 * tmp6 | |
| 285 | 168832264 | tmp8 = tmp7 * tmp7 | |
| 286 | 168832264 | tmp9 = tmp7 * tmp6 | |
| 287 | ans = sqrt(tmp3 * this%n0**6.0_wp - tmp3 * tmp5 * tmp7 + 4.0_wp * tmp3 * tmp5 - & | ||
| 288 | tmp3 * tmp4 * tmp8 + 4.0_wp * tmp3 * tmp4 * tmp7 + tmp3 * tmp8 + 4.0_wp * tmp2 * tmp5 * & | ||
| 289 | tmp6 + 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * tmp9 + tmp1 * tmp5 + 3.0_wp * & | ||
| 290 | tmp1 * tmp4 * tmp7 + 4.0_wp * tmp1 * tmp4 + 6.0_wp * tmp1 * tmp7 + 2.0_wp * & | ||
| 291 | this%R0_pert * tmp4 * tmp6 + 4.0_wp * this%R0_pert * tmp6 + 1.0_wp) / (tmp1 * & | ||
| 292 | 168832264 | tmp4 + tmp1 * tmp7 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp)**(3.0_wp / 2.0_wp) | |
| 293 | 168832264 | end function frenet_serret_curvature | |
| 294 | |||
| 295 | 200 | pure real(wp) function frenet_serret_curvature_tbp(this, zp) result(ans) | |
| 296 | !$acc routine seq | ||
| 297 | implicit none | ||
| 298 | class(FrenetSerretTransform), intent(in) :: this | ||
| 299 | real(wp), intent(in) :: zp | ||
| 300 | 200 | ans = frenet_serret_curvature(this, zp) | |
| 301 | 200 | end function frenet_serret_curvature_tbp | |
| 302 | |||
| 303 | 146699840 | pure real(wp) function frenet_serret_torsion(this, zp) result(ans) | |
| 304 | !$acc routine seq | ||
| 305 | implicit none | ||
| 306 | type(FrenetSerretTransform), intent(in) :: this | ||
| 307 | real(wp), intent(in) :: zp | ||
| 308 | real(wp) :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9 | ||
| 309 | |||
| 310 | 146699840 | tmp1 = this%R0_pert**2.0_wp | |
| 311 | 146699840 | tmp2 = tmp1 * this%R0_pert | |
| 312 | 146699840 | tmp3 = tmp2 * this%R0_pert | |
| 313 | 146699840 | tmp4 = this%n0**2.0_wp | |
| 314 | 146699840 | tmp5 = tmp4 * tmp4 | |
| 315 | 146699840 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 316 | 146699840 | tmp7 = tmp6 * tmp6 | |
| 317 | 146699840 | tmp8 = tmp7 * tmp7 | |
| 318 | 146699840 | tmp9 = tmp7 * tmp6 | |
| 319 | ans = this%R0_pert * this%n0 * (-2.0_wp * tmp1 * tmp5 * tmp6 + tmp1 * tmp4 * tmp9 - & | ||
| 320 | 4.0_wp * tmp1 * tmp4 * tmp6 - tmp1 * tmp9 + this%R0_pert * tmp5 - 4.0_wp * & | ||
| 321 | this%R0_pert * tmp4 * tmp7 + 2.0_wp * this%R0_pert * tmp4 - 2.0_wp * & | ||
| 322 | this%R0_pert * tmp7 + tmp4 * tmp6 - tmp6) / (tmp3 * this%n0**6.0_wp - tmp3 * & | ||
| 323 | tmp5 * tmp7 + 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * tmp8 + 4.0_wp * tmp3 * tmp4 * tmp7 + & | ||
| 324 | tmp3 * tmp8 + 4.0_wp * tmp2 * tmp5 * tmp6 + 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * & | ||
| 325 | tmp2 * tmp9 + tmp1 * tmp5 + 3.0_wp * tmp1 * tmp4 * tmp7 + 4.0_wp * tmp1 * tmp4 + & | ||
| 326 | 6.0_wp * tmp1 * tmp7 + 2.0_wp * this%R0_pert * tmp4 * tmp6 + 4.0_wp * & | ||
| 327 | 146699840 | this%R0_pert * tmp6 + 1.0_wp) | |
| 328 | 146699840 | end function frenet_serret_torsion | |
| 329 | |||
| 330 | 200 | pure real(wp) function frenet_serret_torsion_tbp(this, zp) result(ans) | |
| 331 | !$acc routine seq | ||
| 332 | implicit none | ||
| 333 | class(FrenetSerretTransform), intent(in) :: this | ||
| 334 | real(wp), intent(in) :: zp | ||
| 335 | 200 | ans = frenet_serret_torsion(this, zp) | |
| 336 | 200 | end function frenet_serret_torsion_tbp | |
| 337 | |||
| 338 | ! -- Frenet-Serret frame vectors -------------------------------------- | ||
| 339 | 8700 | pure real(wp) function frenet_serret_tangent(this, i, zp) result(ans) | |
| 340 | !$acc routine seq | ||
| 341 | implicit none | ||
| 342 | type(FrenetSerretTransform), intent(in) :: this | ||
| 343 | integer, intent(in) :: i | ||
| 344 | real(wp), intent(in) :: zp | ||
| 345 | real(wp) :: tmp1, tmp2, tmp3 | ||
| 346 | 11600 | select case (i) | |
| 347 | case (1) | ||
| 348 | 2900 | tmp1 = this%R0_pert**2.0_wp | |
| 349 | 2900 | tmp2 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 350 | 2900 | tmp3 = tmp2 * tmp2 | |
| 351 | ans = (this%R0_pert * this%n0 * cos((6.283185307179586 * zp)) * cos(this%n0 * & | ||
| 352 | (6.283185307179586 * zp)) - (this%R0_pert * tmp2 + 1.0_wp) * & | ||
| 353 | sin((6.283185307179586 * zp))) / sqrt(tmp1 * this%n0**2.0_wp + tmp1 * tmp3 + & | ||
| 354 | 2900 | 2.0_wp * this%R0_pert * tmp2 + 1.0_wp) | |
| 355 | case (2) | ||
| 356 | 2900 | tmp1 = this%R0_pert**2.0_wp | |
| 357 | 2900 | tmp2 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 358 | 2900 | tmp3 = tmp2 * tmp2 | |
| 359 | ans = (this%R0_pert * this%n0 * sin((6.283185307179586 * zp)) * cos(this%n0 * & | ||
| 360 | (6.283185307179586 * zp)) + (this%R0_pert * tmp2 + 1.0_wp) * & | ||
| 361 | cos((6.283185307179586 * zp))) / sqrt(tmp1 * this%n0**2.0_wp + tmp1 * tmp3 + & | ||
| 362 | 2900 | 2.0_wp * this%R0_pert * tmp2 + 1.0_wp) | |
| 363 | case (3) | ||
| 364 | 2900 | tmp1 = this%R0_pert**2.0_wp | |
| 365 | 2900 | tmp2 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 366 | 2900 | tmp3 = tmp2 * tmp2 | |
| 367 | ans = -this%R0_pert * this%n0 * tmp2 / sqrt(tmp1 * this%n0**2.0_wp + tmp1 * & | ||
| 368 |
3/4✓ Branch 2 → 3 taken 2900 times.
✓ Branch 2 → 4 taken 2900 times.
✓ Branch 2 → 5 taken 2900 times.
✗ Branch 2 → 6 not taken.
|
8700 | tmp3 + 2.0_wp * this%R0_pert * tmp2 + 1.0_wp) |
| 369 | end select | ||
| 370 | 8700 | end function frenet_serret_tangent | |
| 371 | |||
| 372 | 1500 | pure real(wp) function frenet_serret_tangent_tbp(this, i, zp) result(ans) | |
| 373 | !$acc routine seq | ||
| 374 | implicit none | ||
| 375 | class(FrenetSerretTransform), intent(in) :: this | ||
| 376 | integer, intent(in) :: i | ||
| 377 | real(wp), intent(in) :: zp | ||
| 378 | 1500 | ans = frenet_serret_tangent(this, i, zp) | |
| 379 | 1500 | end function frenet_serret_tangent_tbp | |
| 380 | |||
| 381 | 10800 | pure real(wp) function frenet_serret_normal(this, i, zp) result(ans) | |
| 382 | !$acc routine seq | ||
| 383 | implicit none | ||
| 384 | type(FrenetSerretTransform), intent(in) :: this | ||
| 385 | integer, intent(in) :: i | ||
| 386 | real(wp), intent(in) :: zp | ||
| 387 | real(wp) :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, tmp12 | ||
| 388 | 14400 | select case (i) | |
| 389 | case (1) | ||
| 390 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 391 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 392 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 393 | 3600 | tmp4 = this%n0**2.0_wp | |
| 394 | 3600 | tmp5 = tmp4 * tmp4 | |
| 395 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 396 | 3600 | tmp7 = cos(this%n0 * (6.283185307179586 * zp)) | |
| 397 | 3600 | tmp8 = cos((6.283185307179586 * zp)) | |
| 398 | 3600 | tmp9 = sin((6.283185307179586 * zp)) | |
| 399 | 3600 | tmp10 = tmp6 * tmp6 | |
| 400 | 3600 | tmp11 = tmp10 * tmp10 | |
| 401 | 3600 | tmp12 = tmp10 * tmp6 | |
| 402 | ans = (-tmp1 * tmp4 * (this%R0_pert * tmp4 * tmp8 + this%R0_pert * this%n0 * & | ||
| 403 | (cos((6.283185307179586 * zp) * (2.0_wp * this%n0 - 1.0_wp)) - & | ||
| 404 | cos((6.283185307179586 * zp) * (2.0_wp * this%n0 + 1.0_wp))) / 4.0_wp + & | ||
| 405 | this%R0_pert * tmp10 * tmp8 - this%n0 * tmp9 * tmp7 + tmp6 * tmp8) * tmp6 - & | ||
| 406 | (this%R0_pert * this%n0 * tmp9 * tmp7 + (this%R0_pert * tmp6 + 1.0_wp) * tmp8) & | ||
| 407 | * (-tmp1 * tmp4 * tmp10 + 2.0_wp * tmp1 * tmp4 + tmp1 * tmp10 + this%R0_pert * & | ||
| 408 | tmp4 * tmp6 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp)) / (sqrt(tmp1 * tmp4 + & | ||
| 409 | tmp1 * tmp10 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp) * sqrt(tmp3 * & | ||
| 410 | this%n0**6.0_wp - tmp3 * tmp5 * tmp10 + 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * & | ||
| 411 | tmp11 + 4.0_wp * tmp3 * tmp4 * tmp10 + tmp3 * tmp11 + 4.0_wp * tmp2 * tmp5 * tmp6 + & | ||
| 412 | 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * tmp12 + tmp1 * tmp5 + 3.0_wp * tmp1 * & | ||
| 413 | tmp4 * tmp10 + 4.0_wp * tmp1 * tmp4 + 6.0_wp * tmp1 * tmp10 + 2.0_wp * & | ||
| 414 | 3600 | this%R0_pert * tmp4 * tmp6 + 4.0_wp * this%R0_pert * tmp6 + 1.0_wp)) | |
| 415 | case (2) | ||
| 416 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 417 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 418 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 419 | 3600 | tmp4 = this%n0**2.0_wp | |
| 420 | 3600 | tmp5 = tmp4 * tmp4 | |
| 421 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 422 | 3600 | tmp7 = cos(this%n0 * (6.283185307179586 * zp)) | |
| 423 | 3600 | tmp8 = cos((6.283185307179586 * zp)) | |
| 424 | 3600 | tmp9 = sin((6.283185307179586 * zp)) | |
| 425 | 3600 | tmp10 = tmp6 * tmp6 | |
| 426 | 3600 | tmp11 = tmp10 * tmp10 | |
| 427 | 3600 | tmp12 = tmp10 * tmp6 | |
| 428 | ans = (-tmp1 * tmp4 * (this%R0_pert * tmp4 * tmp9 - this%R0_pert * this%n0 * & | ||
| 429 | (sin((6.283185307179586 * zp) * (2.0_wp * this%n0 - 1.0_wp)) + & | ||
| 430 | sin((6.283185307179586 * zp) * (2.0_wp * this%n0 + 1.0_wp))) / 4.0_wp + & | ||
| 431 | this%R0_pert * tmp9 * tmp10 + this%n0 * tmp8 * tmp7 + tmp9 * tmp6) * tmp6 + & | ||
| 432 | (this%R0_pert * this%n0 * tmp8 * tmp7 - (this%R0_pert * tmp6 + 1.0_wp) * tmp9) & | ||
| 433 | * (-tmp1 * tmp4 * tmp10 + 2.0_wp * tmp1 * tmp4 + tmp1 * tmp10 + this%R0_pert * & | ||
| 434 | tmp4 * tmp6 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp)) / (sqrt(tmp1 * tmp4 + & | ||
| 435 | tmp1 * tmp10 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp) * sqrt(tmp3 * & | ||
| 436 | this%n0**6.0_wp - tmp3 * tmp5 * tmp10 + 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * & | ||
| 437 | tmp11 + 4.0_wp * tmp3 * tmp4 * tmp10 + tmp3 * tmp11 + 4.0_wp * tmp2 * tmp5 * tmp6 + & | ||
| 438 | 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * tmp12 + tmp1 * tmp5 + 3.0_wp * tmp1 * & | ||
| 439 | tmp4 * tmp10 + 4.0_wp * tmp1 * tmp4 + 6.0_wp * tmp1 * tmp10 + 2.0_wp * & | ||
| 440 | 3600 | this%R0_pert * tmp4 * tmp6 + 4.0_wp * this%R0_pert * tmp6 + 1.0_wp)) | |
| 441 | case (3) | ||
| 442 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 443 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 444 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 445 | 3600 | tmp4 = this%n0**2.0_wp | |
| 446 | 3600 | tmp5 = tmp4 * tmp4 | |
| 447 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 448 | 3600 | tmp7 = tmp6 * tmp6 | |
| 449 | 3600 | tmp8 = tmp7 * tmp7 | |
| 450 | 3600 | tmp9 = tmp7 * tmp6 | |
| 451 | ans = -this%R0_pert * tmp4 * (tmp1 * tmp4 + this%R0_pert * tmp6 + 1.0_wp) * & | ||
| 452 | cos(this%n0 * (6.283185307179586 * zp)) / (sqrt(tmp1 * tmp4 + tmp1 * tmp7 + & | ||
| 453 | 2.0_wp * this%R0_pert * tmp6 + 1.0_wp) * sqrt(tmp3 * this%n0**6.0_wp - tmp3 * & | ||
| 454 | tmp5 * tmp7 + 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * tmp8 + 4.0_wp * tmp3 * tmp4 * & | ||
| 455 | tmp7 + tmp3 * tmp8 + 4.0_wp * tmp2 * tmp5 * tmp6 + 8.0_wp * tmp2 * tmp4 * tmp6 + & | ||
| 456 | 4.0_wp * tmp2 * tmp9 + tmp1 * tmp5 + 3.0_wp * tmp1 * tmp4 * tmp7 + 4.0_wp * tmp1 * & | ||
| 457 | tmp4 + 6.0_wp * tmp1 * tmp7 + 2.0_wp * this%R0_pert * tmp4 * tmp6 + 4.0_wp * & | ||
| 458 |
3/4✓ Branch 2 → 3 taken 3600 times.
✓ Branch 2 → 4 taken 3600 times.
✓ Branch 2 → 5 taken 3600 times.
✗ Branch 2 → 6 not taken.
|
10800 | this%R0_pert * tmp6 + 1.0_wp)) |
| 459 | end select | ||
| 460 | 10800 | end function frenet_serret_normal | |
| 461 | |||
| 462 | 1500 | pure real(wp) function frenet_serret_normal_tbp(this, i, zp) result(ans) | |
| 463 | !$acc routine seq | ||
| 464 | implicit none | ||
| 465 | class(FrenetSerretTransform), intent(in) :: this | ||
| 466 | integer, intent(in) :: i | ||
| 467 | real(wp), intent(in) :: zp | ||
| 468 | 1500 | ans = frenet_serret_normal(this, i, zp) | |
| 469 | 1500 | end function frenet_serret_normal_tbp | |
| 470 | |||
| 471 | 10800 | pure real(wp) function frenet_serret_binormal(this, i, zp) result(ans) | |
| 472 | !$acc routine seq | ||
| 473 | implicit none | ||
| 474 | type(FrenetSerretTransform), intent(in) :: this | ||
| 475 | integer, intent(in) :: i | ||
| 476 | real(wp), intent(in) :: zp | ||
| 477 | real(wp) :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, tmp12 | ||
| 478 | 14400 | select case (i) | |
| 479 | case (1) | ||
| 480 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 481 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 482 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 483 | 3600 | tmp4 = this%n0**2.0_wp | |
| 484 | 3600 | tmp5 = tmp4 * tmp4 | |
| 485 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 486 | 3600 | tmp7 = cos(this%n0 * (6.283185307179586 * zp)) | |
| 487 | 3600 | tmp8 = cos((6.283185307179586 * zp)) | |
| 488 | 3600 | tmp9 = sin((6.283185307179586 * zp)) | |
| 489 | 3600 | tmp10 = tmp6 * tmp6 | |
| 490 | 3600 | tmp11 = tmp10 * tmp10 | |
| 491 | 3600 | tmp12 = tmp10 * tmp6 | |
| 492 | ans = this%R0_pert * this%n0 * (-this%R0_pert * tmp4 * tmp9 + this%R0_pert * & | ||
| 493 | this%n0 * tmp6 * tmp8 * tmp7 - this%R0_pert * tmp9 * tmp10 - this%n0 * tmp8 * & | ||
| 494 | tmp7 - tmp9 * tmp6) / sqrt(tmp3 * this%n0**6.0_wp - tmp3 * tmp5 * tmp10 + & | ||
| 495 | 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * tmp11 + 4.0_wp * tmp3 * tmp4 * tmp10 + tmp3 * & | ||
| 496 | tmp11 + 4.0_wp * tmp2 * tmp5 * tmp6 + 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * & | ||
| 497 | tmp12 + tmp1 * tmp5 + 3.0_wp * tmp1 * tmp4 * tmp10 + 4.0_wp * tmp1 * tmp4 + & | ||
| 498 | 6.0_wp * tmp1 * tmp10 + 2.0_wp * this%R0_pert * tmp4 * tmp6 + 4.0_wp * & | ||
| 499 | 3600 | this%R0_pert * tmp6 + 1.0_wp) | |
| 500 | case (2) | ||
| 501 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 502 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 503 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 504 | 3600 | tmp4 = this%n0**2.0_wp | |
| 505 | 3600 | tmp5 = tmp4 * tmp4 | |
| 506 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 507 | 3600 | tmp7 = cos(this%n0 * (6.283185307179586 * zp)) | |
| 508 | 3600 | tmp8 = cos((6.283185307179586 * zp)) | |
| 509 | 3600 | tmp9 = sin((6.283185307179586 * zp)) | |
| 510 | 3600 | tmp10 = tmp6 * tmp6 | |
| 511 | 3600 | tmp11 = tmp10 * tmp10 | |
| 512 | 3600 | tmp12 = tmp10 * tmp6 | |
| 513 | ans = this%R0_pert * this%n0 * (this%R0_pert * tmp4 * tmp8 + this%R0_pert * & | ||
| 514 | this%n0 * tmp9 * tmp6 * tmp7 + this%R0_pert * tmp10 * tmp8 - this%n0 * tmp9 * & | ||
| 515 | tmp7 + tmp6 * tmp8) / sqrt(tmp3 * this%n0**6.0_wp - tmp3 * tmp5 * tmp10 + & | ||
| 516 | 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * tmp11 + 4.0_wp * tmp3 * tmp4 * tmp10 + tmp3 * & | ||
| 517 | tmp11 + 4.0_wp * tmp2 * tmp5 * tmp6 + 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * & | ||
| 518 | tmp12 + tmp1 * tmp5 + 3.0_wp * tmp1 * tmp4 * tmp10 + 4.0_wp * tmp1 * tmp4 + & | ||
| 519 | 6.0_wp * tmp1 * tmp10 + 2.0_wp * this%R0_pert * tmp4 * tmp6 + 4.0_wp * & | ||
| 520 | 3600 | this%R0_pert * tmp6 + 1.0_wp) | |
| 521 | case (3) | ||
| 522 | 3600 | tmp1 = this%R0_pert**2.0_wp | |
| 523 | 3600 | tmp2 = tmp1 * this%R0_pert | |
| 524 | 3600 | tmp3 = tmp2 * this%R0_pert | |
| 525 | 3600 | tmp4 = this%n0**2.0_wp | |
| 526 | 3600 | tmp5 = tmp4 * tmp4 | |
| 527 | 3600 | tmp6 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 528 | 3600 | tmp7 = tmp6 * tmp6 | |
| 529 | 3600 | tmp8 = tmp7 * tmp7 | |
| 530 | 3600 | tmp9 = tmp7 * tmp6 | |
| 531 | ans = (-tmp1 * tmp4 * tmp7 + 2.0_wp * tmp1 * tmp4 + tmp1 * tmp7 + this%R0_pert * & | ||
| 532 | tmp4 * tmp6 + 2.0_wp * this%R0_pert * tmp6 + 1.0_wp) / sqrt(tmp3 * & | ||
| 533 | this%n0**6.0_wp - tmp3 * tmp5 * tmp7 + 4.0_wp * tmp3 * tmp5 - tmp3 * tmp4 * & | ||
| 534 | tmp8 + 4.0_wp * tmp3 * tmp4 * tmp7 + tmp3 * tmp8 + 4.0_wp * tmp2 * tmp5 * tmp6 + & | ||
| 535 | 8.0_wp * tmp2 * tmp4 * tmp6 + 4.0_wp * tmp2 * tmp9 + tmp1 * tmp5 + 3.0_wp * tmp1 * & | ||
| 536 | tmp4 * tmp7 + 4.0_wp * tmp1 * tmp4 + 6.0_wp * tmp1 * tmp7 + 2.0_wp * & | ||
| 537 |
3/4✓ Branch 2 → 3 taken 3600 times.
✓ Branch 2 → 4 taken 3600 times.
✓ Branch 2 → 5 taken 3600 times.
✗ Branch 2 → 6 not taken.
|
10800 | this%R0_pert * tmp4 * tmp6 + 4.0_wp * this%R0_pert * tmp6 + 1.0_wp) |
| 538 | end select | ||
| 539 | 10800 | end function frenet_serret_binormal | |
| 540 | |||
| 541 | 1500 | pure real(wp) function frenet_serret_binormal_tbp(this, i, zp) result(ans) | |
| 542 | !$acc routine seq | ||
| 543 | implicit none | ||
| 544 | class(FrenetSerretTransform), intent(in) :: this | ||
| 545 | integer, intent(in) :: i | ||
| 546 | real(wp), intent(in) :: zp | ||
| 547 | 1500 | ans = frenet_serret_binormal(this, i, zp) | |
| 548 | 1500 | end function frenet_serret_binormal_tbp | |
| 549 | |||
| 550 | 168832164 | pure real(wp) function frenet_serret_derivative_magnitude(this, zp) result(ans) | |
| 551 | !$acc routine seq | ||
| 552 | implicit none | ||
| 553 | type(FrenetSerretTransform), intent(in) :: this | ||
| 554 | real(wp), intent(in) :: zp | ||
| 555 | real(wp) :: tmp1, tmp2, tmp3 | ||
| 556 | |||
| 557 | 168832164 | tmp1 = this%R0_pert**2.0_wp | |
| 558 | 168832164 | tmp2 = sin(this%n0 * (6.283185307179586 * zp)) | |
| 559 | 168832164 | tmp3 = tmp2 * tmp2 | |
| 560 | ans = sqrt(tmp1 * this%n0**2.0_wp + tmp1 * tmp3 + 2.0_wp * this%R0_pert * tmp2 + & | ||
| 561 | 168832164 | 1.0_wp) | |
| 562 | 168832164 | ans = ans * 6.283185307179586 ! Apply scaling factor for z-coordinate | |
| 563 | 168832164 | end function frenet_serret_derivative_magnitude | |
| 564 | |||
| 565 | 100 | pure real(wp) function frenet_serret_derivative_magnitude_tbp(this, zp) result(ans) | |
| 566 | !$acc routine seq | ||
| 567 | implicit none | ||
| 568 | class(FrenetSerretTransform), intent(in) :: this | ||
| 569 | real(wp), intent(in) :: zp | ||
| 570 | 100 | ans = frenet_serret_derivative_magnitude(this, zp) | |
| 571 | 100 | end function frenet_serret_derivative_magnitude_tbp | |
| 572 | |||
| 573 | 22132424 | pure real(wp) function frenet_serret_det2D(this, xp, yp, zp) result(ans) | |
| 574 | !$acc routine seq | ||
| 575 | implicit none | ||
| 576 | type(FrenetSerretTransform), intent(in) :: this | ||
| 577 | real(wp), intent(in) :: xp, yp, zp | ||
| 578 | |||
| 579 | ans = this%R1**2.0_wp * xp * (this%R1_pert * xp * cos(this%n1 * (6.283185307179586 * & | ||
| 580 | yp)) + 1.0_wp) * (2.0_wp * this%R1_pert * xp * cos(this%n1 * (6.283185307179586 * & | ||
| 581 | 22132424 | yp)) + 1.0_wp) | |
| 582 | 22132424 | ans = ans * 6.283185307179586 ! Apply scaling factor for y-coordinate | |
| 583 | 22132424 | end function frenet_serret_det2D | |
| 584 | |||
| 585 | ✗ | pure real(wp) function frenet_serret_inverse_transform(this, i, x, y, z) result(ans) | |
| 586 | !$acc routine seq | ||
| 587 | implicit none | ||
| 588 | type(FrenetSerretTransform), intent(in) :: this | ||
| 589 | integer, intent(in) :: i | ||
| 590 | real(wp), intent(in) :: x, y, z | ||
| 591 | character(len=256) :: error_msg | ||
| 592 | ! TODO error handling on GPU | ||
| 593 | ! write (error_msg, '(A)') 'ERROR in FrenetSerretTransform::inverse_transform: inverse transform is not implemented.' | ||
| 594 | ! error stop error_msg | ||
| 595 | ans = 0._wp | ||
| 596 | ✗ | end function frenet_serret_inverse_transform | |
| 597 | |||
| 598 | ! -- Decomposition functions R, C, S --------------------------------- | ||
| 599 | 168834164 | pure real(wp) function frenet_serret_R(this, xp, yp, zp) result(ans) | |
| 600 | !$acc routine seq | ||
| 601 | implicit none | ||
| 602 | type(FrenetSerretTransform), intent(in) :: this | ||
| 603 | real(wp), intent(in) :: xp, yp, zp | ||
| 604 | ans = this%R1 * xp * (this%R1_pert * xp * cos(this%n1 * (6.283185307179586 * yp)) + & | ||
| 605 | 168834164 | 1.0_wp) | |
| 606 | 168834164 | end function frenet_serret_R | |
| 607 | |||
| 608 | pure real(wp) function frenet_serret_C(this, xp, yp, zp) result(ans) | ||
| 609 | !$acc routine seq | ||
| 610 | implicit none | ||
| 611 | type(FrenetSerretTransform), intent(in) :: this | ||
| 612 | real(wp), intent(in) :: xp, yp, zp | ||
| 613 | 146701740 | ans = cos((6.283185307179586 * yp)) | |
| 614 | end function frenet_serret_C | ||
| 615 | |||
| 616 | pure real(wp) function frenet_serret_S(this, xp, yp, zp) result(ans) | ||
| 617 | !$acc routine seq | ||
| 618 | implicit none | ||
| 619 | type(FrenetSerretTransform), intent(in) :: this | ||
| 620 | real(wp), intent(in) :: xp, yp, zp | ||
| 621 | 168834164 | ans = sin((6.283185307179586 * yp)) | |
| 622 | end function frenet_serret_S | ||
| 623 | |||
| 624 | 146699640 | pure real(wp) function frenet_serret_Rx(this, xp, yp, zp) result(ans) | |
| 625 | !$acc routine seq | ||
| 626 | implicit none | ||
| 627 | type(FrenetSerretTransform), intent(in) :: this | ||
| 628 | real(wp), intent(in) :: xp, yp, zp | ||
| 629 | ans = this%R1 * (2.0_wp * this%R1_pert * xp * cos(this%n1 * (6.283185307179586 * yp)) & | ||
| 630 | 146699640 | + 1.0_wp) | |
| 631 | 146699640 | end function frenet_serret_Rx | |
| 632 | |||
| 633 | 146699640 | pure real(wp) function frenet_serret_Ry(this, xp, yp, zp) result(ans) | |
| 634 | !$acc routine seq | ||
| 635 | implicit none | ||
| 636 | type(FrenetSerretTransform), intent(in) :: this | ||
| 637 | real(wp), intent(in) :: xp, yp, zp | ||
| 638 | ans = -this%R1 * this%R1_pert * this%n1 * xp**2.0_wp * sin(this%n1 * & | ||
| 639 | 146699640 | (6.283185307179586 * yp)) | |
| 640 | 146699640 | ans = ans * 6.283185307179586 ! Apply scaling factor for y-coordinate | |
| 641 | 146699640 | end function frenet_serret_Ry | |
| 642 | |||
| 643 | pure real(wp) function frenet_serret_Cy(this, xp, yp, zp) result(ans) | ||
| 644 | !$acc routine seq | ||
| 645 | implicit none | ||
| 646 | type(FrenetSerretTransform), intent(in) :: this | ||
| 647 | real(wp), intent(in) :: xp, yp, zp | ||
| 648 | 146699640 | ans = -sin((6.283185307179586 * yp)) | |
| 649 | 146699640 | ans = ans * 6.283185307179586 ! Apply scaling factor for y-coordinate | |
| 650 | end function frenet_serret_Cy | ||
| 651 | |||
| 652 | pure real(wp) function frenet_serret_Sy(this, xp, yp, zp) result(ans) | ||
| 653 | !$acc routine seq | ||
| 654 | implicit none | ||
| 655 | type(FrenetSerretTransform), intent(in) :: this | ||
| 656 | real(wp), intent(in) :: xp, yp, zp | ||
| 657 | ans = cos((6.283185307179586 * yp)) | ||
| 658 |
4/5✓ Branch 2 → 3 taken 88017384 times.
✓ Branch 2 → 4 taken 2400 times.
✓ Branch 2 → 5 taken 2400 times.
✗ Branch 2 → 6 not taken.
✓ Branch 2 → 7 taken 58677456 times.
|
146699640 | ans = ans * 6.283185307179586 ! Apply scaling factor for y-coordinate |
| 659 | end function frenet_serret_Sy | ||
| 660 | |||
| 661 | ! -- Base curve and derivative ---------------------------------------- | ||
| 662 | 2100 | pure real(wp) function frenet_serret_curve(this, i, zp) result(ans) | |
| 663 | !$acc routine seq | ||
| 664 | implicit none | ||
| 665 | type(FrenetSerretTransform), intent(in) :: this | ||
| 666 | integer, intent(in) :: i | ||
| 667 | real(wp), intent(in) :: zp | ||
| 668 | 2800 | select case (i) | |
| 669 | case (1) | ||
| 670 | ans = (this%R0_pert * sin(this%n0 * (6.283185307179586 * zp)) + 1.0_wp) * & | ||
| 671 | 700 | cos((6.283185307179586 * zp)) | |
| 672 | case (2) | ||
| 673 | ans = (this%R0_pert * sin(this%n0 * (6.283185307179586 * zp)) + 1.0_wp) * & | ||
| 674 | 700 | sin((6.283185307179586 * zp)) | |
| 675 | case (3) | ||
| 676 |
3/4✓ Branch 2 → 3 taken 700 times.
✓ Branch 2 → 4 taken 700 times.
✓ Branch 2 → 5 taken 700 times.
✗ Branch 2 → 6 not taken.
|
2100 | ans = this%R0_pert * cos(this%n0 * (6.283185307179586 * zp)) |
| 677 | end select | ||
| 678 | 2100 | end function frenet_serret_curve | |
| 679 | |||
| 680 | 2400 | pure real(wp) function frenet_serret_curve_derivative(this, i, zp) result(ans) | |
| 681 | !$acc routine seq | ||
| 682 | implicit none | ||
| 683 | type(FrenetSerretTransform), intent(in) :: this | ||
| 684 | integer, intent(in) :: i | ||
| 685 | real(wp), intent(in) :: zp | ||
| 686 | 3200 | select case (i) | |
| 687 | case (1) | ||
| 688 | ans = this%R0_pert * this%n0 * cos((6.283185307179586 * zp)) * cos(this%n0 * & | ||
| 689 | (6.283185307179586 * zp)) - (this%R0_pert * sin(this%n0 * & | ||
| 690 | 800 | (6.283185307179586 * zp)) + 1.0_wp) * sin((6.283185307179586 * zp)) | |
| 691 | case (2) | ||
| 692 | ans = this%R0_pert * this%n0 * sin((6.283185307179586 * zp)) * cos(this%n0 * & | ||
| 693 | (6.283185307179586 * zp)) + (this%R0_pert * sin(this%n0 * & | ||
| 694 | 800 | (6.283185307179586 * zp)) + 1.0_wp) * cos((6.283185307179586 * zp)) | |
| 695 | case (3) | ||
| 696 |
3/4✓ Branch 2 → 3 taken 800 times.
✓ Branch 2 → 4 taken 800 times.
✓ Branch 2 → 5 taken 800 times.
✗ Branch 2 → 6 not taken.
|
2400 | ans = -this%R0_pert * this%n0 * sin(this%n0 * (6.283185307179586 * zp)) |
| 697 | end select | ||
| 698 | 2400 | ans = ans * 6.283185307179586 ! Apply scaling factor for z-coordinate | |
| 699 | 2400 | end function frenet_serret_curve_derivative | |
| 700 | |||
| 701 | ✗ | end module m_coord_transform_frenet_serret | |
| 702 | |||
| 703 |