bessel/m_bessel.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module for Bessel functions and their zeros | ||
| 2 | module m_bessel | ||
| 3 | use m_bessel_roots_data, only: MAX_ORDER, MAX_NR_ROOTS | ||
| 4 | use m_common, only: wp, pi | ||
| 5 | |||
| 6 | private | ||
| 7 | public :: bessel_jn_zero, bessel_jn_derivative, bessel_jn_prime_zero, MAX_ORDER, MAX_NR_ROOTS | ||
| 8 | |||
| 9 | public :: disc_laplace_eigenvalues, disc_curlcurl_eigenvalues | ||
| 10 | |||
| 11 | public :: init_manufactured_solution_bessel | ||
| 12 | public :: toroidal_zeroform | ||
| 13 | public :: toroidal_oneform_x, toroidal_oneform_y, toroidal_oneform_z | ||
| 14 | public :: toroidal_twoform_x, toroidal_twoform_y, toroidal_twoform_z | ||
| 15 | public :: toroidal_threeform | ||
| 16 | |||
| 17 | ! Parameters for the manufactured toroidal solution | ||
| 18 | integer :: param_bessel_n = -1 ! Order of the Bessel function | ||
| 19 | integer :: param_bessel_m = -1 ! Index of the zero | ||
| 20 | real(wp) :: lambda_bessel_mn = -1._wp ! m-th zero of the Bessel function of order n | ||
| 21 | real(wp) :: coords_major_radius = 1.0_wp ! Major radius of the torus | ||
| 22 | real(wp) :: coords_minor_radius = 0.2_wp ! Minor radius of the torus | ||
| 23 | |||
| 24 | real(wp) :: lambda_target_ | ||
| 25 | contains | ||
| 26 | !> @brief Compute the m-th zero of the Bessel function of the first kind of order n | ||
| 27 | !> | ||
| 28 | !> @param[in] n The order of the Bessel function | ||
| 29 | !> @param[in] m The index of the zero | ||
| 30 | !> | ||
| 31 | !> @return The m-th zero of the Bessel function of the first kind of order n | ||
| 32 | 984 | pure real(wp) function bessel_jn_zero(n, m) | |
| 33 | use m_bessel_roots_data | ||
| 34 | |||
| 35 | implicit none | ||
| 36 | integer, intent(in) :: n, m | ||
| 37 | character(len=256) :: error_msg | ||
| 38 | |||
| 39 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 11 taken 984 times.
|
984 | if (m < 1 .or. m > MAX_NR_ROOTS) then |
| 40 | ✗ | write (error_msg, '(A,I0,A,I0,A)') "ERROR in bessel_jn_zero: root index m = ", m, & | |
| 41 | ✗ | " exceeds maximum number of roots (MAX_NR_ROOTS = ", MAX_NR_ROOTS, ")" | |
| 42 | ✗ | error stop error_msg | |
| 43 | end if | ||
| 44 |
1/2✗ Branch 11 → 12 not taken.
✓ Branch 11 → 19 taken 984 times.
|
984 | if (n < 0 .or. n > MAX_ORDER) then |
| 45 | ✗ | write (error_msg, '(A,I0,A,I0)') "ERROR in bessel_jn_zero: order n = ", n, & | |
| 46 | ✗ | " is out of range [0, MAX_ORDER]. MAX_ORDER = ", MAX_ORDER | |
| 47 | ✗ | error stop error_msg | |
| 48 | end if | ||
| 49 | |||
| 50 | 64 | select case (n) | |
| 51 | case (0) | ||
| 52 | 64 | bessel_jn_zero = bessel0_roots(m) | |
| 53 | case (1) | ||
| 54 | 116 | bessel_jn_zero = bessel1_roots(m) | |
| 55 | case (2) | ||
| 56 | 102 | bessel_jn_zero = bessel2_roots(m) | |
| 57 | case (3) | ||
| 58 | 102 | bessel_jn_zero = bessel3_roots(m) | |
| 59 | case (4) | ||
| 60 | 100 | bessel_jn_zero = bessel4_roots(m) | |
| 61 | case (5) | ||
| 62 | 100 | bessel_jn_zero = bessel5_roots(m) | |
| 63 | case (6) | ||
| 64 | 80 | bessel_jn_zero = bessel6_roots(m) | |
| 65 | case (7) | ||
| 66 | 80 | bessel_jn_zero = bessel7_roots(m) | |
| 67 | case (8) | ||
| 68 | 80 | bessel_jn_zero = bessel8_roots(m) | |
| 69 | case (9) | ||
| 70 | 80 | bessel_jn_zero = bessel9_roots(m) | |
| 71 | case (10) | ||
| 72 |
11/11✓ Branch 19 → 20 taken 64 times.
✓ Branch 19 → 21 taken 116 times.
✓ Branch 19 → 22 taken 102 times.
✓ Branch 19 → 23 taken 102 times.
✓ Branch 19 → 24 taken 100 times.
✓ Branch 19 → 25 taken 100 times.
✓ Branch 19 → 26 taken 80 times.
✓ Branch 19 → 27 taken 80 times.
✓ Branch 19 → 28 taken 80 times.
✓ Branch 19 → 29 taken 80 times.
✓ Branch 19 → 30 taken 80 times.
|
984 | bessel_jn_zero = bessel10_roots(m) |
| 73 | end select | ||
| 74 | 984 | end function | |
| 75 | |||
| 76 | !> @brief Compute the m-th zero of the derivative of the Bessel function of the first kind of order n | ||
| 77 | !> | ||
| 78 | !> @param[in] n The order of the Bessel function | ||
| 79 | !> @param[in] m The index of the zero | ||
| 80 | !> | ||
| 81 | !> @return The m-th zero of the derivative of the Bessel function of the first kind of order n | ||
| 82 | 421 | pure real(wp) function bessel_jn_prime_zero(n, m) | |
| 83 | use m_bessel_prime_roots_data | ||
| 84 | |||
| 85 | implicit none | ||
| 86 | integer, intent(in) :: n, m | ||
| 87 | character(len=256) :: error_msg | ||
| 88 | |||
| 89 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 11 taken 421 times.
|
421 | if (m < 1 .or. m > MAX_NR_ROOTS) then |
| 90 | ✗ | write (error_msg, '(A,I0,A,I0,A)') "ERROR in bessel_jn_prime_zero: root index m = ", m, & | |
| 91 | ✗ | " exceeds maximum number of roots (MAX_NR_ROOTS = ", MAX_NR_ROOTS, ")" | |
| 92 | ✗ | error stop error_msg | |
| 93 | end if | ||
| 94 |
1/2✗ Branch 11 → 12 not taken.
✓ Branch 11 → 19 taken 421 times.
|
421 | if (n < 0 .or. n > MAX_ORDER) then |
| 95 | ✗ | write (error_msg, '(A,I0,A,I0)') "ERROR in bessel_jn_prime_zero: order n = ", n, & | |
| 96 | ✗ | " is out of range [0, MAX_ORDER]. MAX_ORDER = ", MAX_ORDER | |
| 97 | ✗ | error stop error_msg | |
| 98 | end if | ||
| 99 | |||
| 100 | 21 | select case (n) | |
| 101 | case (0) | ||
| 102 | 21 | bessel_jn_prime_zero = bessel_prime0_roots(m) | |
| 103 | case (1) | ||
| 104 | 40 | bessel_jn_prime_zero = bessel_prime1_roots(m) | |
| 105 | case (2) | ||
| 106 | 40 | bessel_jn_prime_zero = bessel_prime2_roots(m) | |
| 107 | case (3) | ||
| 108 | 40 | bessel_jn_prime_zero = bessel_prime3_roots(m) | |
| 109 | case (4) | ||
| 110 | 40 | bessel_jn_prime_zero = bessel_prime4_roots(m) | |
| 111 | case (5) | ||
| 112 | 40 | bessel_jn_prime_zero = bessel_prime5_roots(m) | |
| 113 | case (6) | ||
| 114 | 40 | bessel_jn_prime_zero = bessel_prime6_roots(m) | |
| 115 | case (7) | ||
| 116 | 40 | bessel_jn_prime_zero = bessel_prime7_roots(m) | |
| 117 | case (8) | ||
| 118 | 40 | bessel_jn_prime_zero = bessel_prime8_roots(m) | |
| 119 | case (9) | ||
| 120 | 40 | bessel_jn_prime_zero = bessel_prime9_roots(m) | |
| 121 | case (10) | ||
| 122 |
11/11✓ Branch 19 → 20 taken 21 times.
✓ Branch 19 → 21 taken 40 times.
✓ Branch 19 → 22 taken 40 times.
✓ Branch 19 → 23 taken 40 times.
✓ Branch 19 → 24 taken 40 times.
✓ Branch 19 → 25 taken 40 times.
✓ Branch 19 → 26 taken 40 times.
✓ Branch 19 → 27 taken 40 times.
✓ Branch 19 → 28 taken 40 times.
✓ Branch 19 → 29 taken 40 times.
✓ Branch 19 → 30 taken 40 times.
|
421 | bessel_jn_prime_zero = bessel_prime10_roots(m) |
| 123 | end select | ||
| 124 | 421 | end function | |
| 125 | |||
| 126 | !> @brief Compute the m-th zero of the Bessel function of the first kind of order n or its derivative | ||
| 127 | !> | ||
| 128 | !> @param[in] n The order of the Bessel function | ||
| 129 | !> @param[in] m The index of the zero | ||
| 130 | !> @param[in] prime Logical flag indicating whether to compute the zero of the derivative (true) or the function itself (false) | ||
| 131 | !> | ||
| 132 | !> @return The m-th zero of the Bessel function of the first kind of order n or its derivative | ||
| 133 | 1263 | pure real(wp) function bessel_jn_ANY_zero(n, m, prime) | |
| 134 | |||
| 135 | implicit none | ||
| 136 | integer, intent(in) :: n, m | ||
| 137 | logical, intent(in) :: prime | ||
| 138 | |||
| 139 |
2/2✓ Branch 2 → 3 taken 421 times.
✓ Branch 2 → 4 taken 842 times.
|
1263 | if (prime) then |
| 140 | 421 | bessel_jn_ANY_zero = bessel_jn_prime_zero(n, m) | |
| 141 | else | ||
| 142 | 842 | bessel_jn_ANY_zero = bessel_jn_zero(n, m) | |
| 143 | end if | ||
| 144 | 1263 | end function | |
| 145 | |||
| 146 | !> @brief Compute the derivative of the Bessel function of the first kind of order n at point x | ||
| 147 | !> | ||
| 148 | !> @param[in] n The order of the Bessel function | ||
| 149 | !> @param[in] x The point at which to evaluate the derivative | ||
| 150 | !> | ||
| 151 | !> @return The value of the derivative | ||
| 152 | 52919796 | pure real(wp) function bessel_jn_derivative(n, x) result(ans) | |
| 153 | implicit none | ||
| 154 | |||
| 155 | integer, intent(in) :: n | ||
| 156 | real(wp), intent(in) :: x | ||
| 157 | |||
| 158 |
1/2✓ Branch 2 → 3 taken 52919796 times.
✗ Branch 2 → 4 not taken.
|
52919796 | if (n > 0) then |
| 159 | 52919796 | ans = (bessel_jn(n - 1, x) - bessel_jn(n + 1, x)) / 2 | |
| 160 | else ! if (n == 0) | ||
| 161 | ✗ | ans = -bessel_jn(1, x) | |
| 162 | end if | ||
| 163 | 52919796 | end function | |
| 164 | |||
| 165 | !> @brief Compute the m-th eigenvalues of the Laplace operator on the unit disk | ||
| 166 | !> | ||
| 167 | !> @param[in] lambda_r The real part of the eigenvalues | ||
| 168 | !> @param[in] n The number of eigenvalues to compute | ||
| 169 | !> @param[in] lambda_target _(optional)_ Target eigenvalue (default: 0) | ||
| 170 | 2 | subroutine disc_laplace_eigenvalues(lambda_r, n, lambda_target) | |
| 171 | implicit none | ||
| 172 | real(wp), allocatable, intent(out) :: lambda_r(:) | ||
| 173 | integer, intent(in) :: n | ||
| 174 | real(wp), intent(in), optional :: lambda_target | ||
| 175 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 2 times.
|
2 | call disc_ANY_eigenvalues(lambda_r, n, .false., lambda_target) |
| 176 | 2 | end subroutine | |
| 177 | |||
| 178 | !> @brief Compute the m-th eigenvalues of the curlcurl operator on the unit disk | ||
| 179 | !> | ||
| 180 | !> @param[in] lambda_r The real part of the eigenvalues | ||
| 181 | !> @param[in] n The number of eigenvalues to compute | ||
| 182 | !> @param[in] lambda_target _(optional)_ Target eigenvalue (default: 0) | ||
| 183 | 1 | subroutine disc_curlcurl_eigenvalues(lambda_r, n, lambda_target) | |
| 184 | implicit none | ||
| 185 | real(wp), allocatable, intent(out) :: lambda_r(:) | ||
| 186 | integer, intent(in) :: n | ||
| 187 | real(wp), intent(in), optional :: lambda_target | ||
| 188 | |||
| 189 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 1 time.
|
1 | call disc_ANY_eigenvalues(lambda_r, n, .true., lambda_target) |
| 190 | 1 | end subroutine | |
| 191 | |||
| 192 | !> @brief Compute the eigenvalues of the Laplace/curlcurl operator on the unit disk | ||
| 193 | !> | ||
| 194 | !> @param[in] lambda_r The real part of the eigenvalues | ||
| 195 | !> @param[in] n The number of eigenvalues to compute | ||
| 196 | !> @param[in] prime Logical flag indicating whether to use zeros of the derivative (true) or the function itself (false) | ||
| 197 | !> @param[in] lambda_target _(optional)_ Target eigenvalue (default: 0) | ||
| 198 | 3 | subroutine disc_ANY_eigenvalues(lambda_r, n, prime, lambda_target) | |
| 199 | use m_bessel_roots_data, only: MAX_NR_ROOTS, MAX_ORDER | ||
| 200 | use m_common, only: sort | ||
| 201 | implicit none | ||
| 202 | |||
| 203 | real(wp), allocatable, intent(out) :: lambda_r(:) | ||
| 204 | integer, intent(in) :: n | ||
| 205 | logical, intent(in) :: prime | ||
| 206 | real(wp), intent(in), optional :: lambda_target | ||
| 207 | |||
| 208 | real(wp) :: all_eigenvalues(2 * MAX_NR_ROOTS * (MAX_ORDER + 1)) | ||
| 209 | integer :: order, m, count, i | ||
| 210 | |||
| 211 | ! Generate all eigenvalues with proper multiplicities | ||
| 212 | 3 | count = 0 | |
| 213 |
2/2✓ Branch 3 → 4 taken 33 times.
✓ Branch 3 → 16 taken 3 times.
|
36 | do order = 0, MAX_ORDER |
| 214 |
2/2✓ Branch 5 → 6 taken 660 times.
✓ Branch 5 → 14 taken 33 times.
|
696 | do m = 1, MAX_NR_ROOTS |
| 215 |
2/2✓ Branch 6 → 7 taken 60 times.
✓ Branch 6 → 9 taken 600 times.
|
693 | if (order == 0) then |
| 216 | ! n=0: multiplicity 1 | ||
| 217 | 60 | count = count + 1 | |
| 218 |
1/2✓ Branch 7 → 8 taken 60 times.
✗ Branch 7 → 13 not taken.
|
60 | if (count <= size(all_eigenvalues)) then |
| 219 | 60 | all_eigenvalues(count) = bessel_jn_ANY_zero(order, m, prime)**2 | |
| 220 | end if | ||
| 221 | else | ||
| 222 | ! n>0: multiplicity 2 | ||
| 223 | 600 | count = count + 1 | |
| 224 |
1/2✓ Branch 9 → 10 taken 600 times.
✗ Branch 9 → 11 not taken.
|
600 | if (count <= size(all_eigenvalues)) then |
| 225 | 600 | all_eigenvalues(count) = bessel_jn_ANY_zero(order, m, prime)**2 | |
| 226 | end if | ||
| 227 | 600 | count = count + 1 | |
| 228 |
1/2✓ Branch 11 → 12 taken 600 times.
✗ Branch 11 → 13 not taken.
|
600 | if (count <= size(all_eigenvalues)) then |
| 229 | 600 | all_eigenvalues(count) = bessel_jn_ANY_zero(order, m, prime)**2 | |
| 230 | end if | ||
| 231 | end if | ||
| 232 | end do | ||
| 233 | end do | ||
| 234 | |||
| 235 | ! Sort eigenvalues by magnitude using the generic sort interface | ||
| 236 |
2/2✓ Branch 17 → 18 taken 1 time.
✓ Branch 17 → 20 taken 2 times.
|
3 | if (present(lambda_target)) then |
| 237 | 1 | lambda_target_ = lambda_target | |
| 238 | 1 | call sort(all_eigenvalues(1:count), sort_fun=sort_fun) | |
| 239 | else | ||
| 240 | 2 | call sort(all_eigenvalues(1:count)) | |
| 241 | end if | ||
| 242 | |||
| 243 |
1/2✗ Branch 22 → 23 not taken.
✓ Branch 22 → 33 taken 3 times.
|
3 | if (n > count) then |
| 244 | ✗ | write (*, '(A,I0,A,I0)') "WARNING in disc_ANY_eigenvalues: requested number of eigenvalues n = ", n, & | |
| 245 | ✗ | " exceeds the number of available eigenvalues (", count, "). Returning only ", count, " eigenvalues." | |
| 246 | end if | ||
| 247 | |||
| 248 | 3 | count = min(n, count) | |
| 249 | |||
| 250 | ! Extract the first n eigenvalues | ||
| 251 |
7/14✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 3 times.
✓ Branch 35 → 34 taken 3 times.
✗ Branch 35 → 36 not taken.
✓ Branch 36 → 37 taken 3 times.
✗ Branch 36 → 38 not taken.
✓ Branch 38 → 39 taken 3 times.
✗ Branch 38 → 40 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 3 times.
✗ Branch 42 → 43 not taken.
✓ Branch 42 → 44 taken 3 times.
✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 3 times.
|
9 | allocate (lambda_r(count)) |
| 252 |
2/2✓ Branch 47 → 48 taken 25 times.
✓ Branch 47 → 49 taken 3 times.
|
28 | do i = 1, count |
| 253 | 28 | lambda_r(i) = all_eigenvalues(i) | |
| 254 | end do | ||
| 255 | |||
| 256 |
1/2✗ Branch 50 → 51 not taken.
✓ Branch 50 → 55 taken 3 times.
|
3 | if (lambda_r(count) > bessel_jn_ANY_zero(0, MAX_NR_ROOTS, prime)**2) then |
| 257 | write (*, '(A)') "WARNING in disc_ANY_eigenvalues: the largest sorted eigenvalue exceeds the square of the largest"// & | ||
| 258 | ✗ | & " stored zeroth-order Bessel root. This means that some of the largest computed eigenvalues may be missing." | |
| 259 | end if | ||
| 260 | contains | ||
| 261 | 5226 | pure real(wp) function sort_fun(x) | |
| 262 | real(wp), intent(in) :: x | ||
| 263 | 5226 | sort_fun = abs(x - lambda_target_) | |
| 264 | 5226 | end function sort_fun | |
| 265 | end subroutine | ||
| 266 | |||
| 267 | !> @brief Initialize the parameters for the manufactured solution involving Bessel functions | ||
| 268 | !> @param[in] n The order of the Bessel function | ||
| 269 | !> @param[in] m The index of the zero | ||
| 270 | !> @param[in] major_radius The major radius of the torus | ||
| 271 | !> @param[in] minor_radius The minor radius of the torus | ||
| 272 | 13 | subroutine init_manufactured_solution_bessel(n, m, major_radius, minor_radius) | |
| 273 | use m_bessel_roots_data, only: MAX_ORDER, MAX_NR_ROOTS | ||
| 274 | implicit none | ||
| 275 | |||
| 276 | integer, intent(in) :: n, m | ||
| 277 | real(wp), intent(in) :: major_radius, minor_radius | ||
| 278 | |||
| 279 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 10 taken 13 times.
|
13 | if (n < 0 .or. n > MAX_ORDER) then |
| 280 | ✗ | write (*, '(A,I0,A,I0)') "ERROR in init_manufactured_solution_bessel: order n = ", n, & | |
| 281 | ✗ | " is out of range [0, MAX_ORDER]. MAX_ORDER = ", MAX_ORDER | |
| 282 | ✗ | error stop 1 | |
| 283 | end if | ||
| 284 |
1/2✗ Branch 10 → 11 not taken.
✓ Branch 10 → 19 taken 13 times.
|
13 | if (m < 1 .or. m > MAX_NR_ROOTS) then |
| 285 | ✗ | write (*, '(A,I0,A,I0)') "ERROR in init_manufactured_solution_bessel: root index m = ", m, & | |
| 286 | ✗ | " exceeds maximum number of roots (MAX_NR_ROOTS = ", MAX_NR_ROOTS, ")" | |
| 287 | ✗ | error stop 1 | |
| 288 | end if | ||
| 289 | |||
| 290 | 13 | param_bessel_n = n | |
| 291 | 13 | param_bessel_m = m | |
| 292 | 13 | lambda_bessel_mn = bessel_jn_zero(n, m) | |
| 293 | 13 | coords_major_radius = major_radius | |
| 294 | 13 | coords_minor_radius = minor_radius | |
| 295 | 13 | end subroutine | |
| 296 | |||
| 297 | !> @brief Compute the part of the manufactured solution involving Bessel functions | ||
| 298 | !> @param[in] xp The first logical coordinate | ||
| 299 | !> @param[in] yp The second logical coordinate | ||
| 300 | !> | ||
| 301 | !> @return The value of the Bessel part of the solution | ||
| 302 | 186754356 | pure real(wp) function bessel_part(xp, yp) | |
| 303 | real(wp), intent(in) :: xp, yp | ||
| 304 | |||
| 305 | bessel_part = bessel_jn(param_bessel_n, lambda_bessel_mn * xp) * & | ||
| 306 | 186754356 | cos(2 * pi * param_bessel_n * yp) | |
| 307 | 186754356 | end function bessel_part | |
| 308 | |||
| 309 | !> @brief Compute the derivative of the Bessel part with respect to x | ||
| 310 | !> @param[in] xp The first logical coordinate | ||
| 311 | !> @param[in] yp The second logical coordinate | ||
| 312 | !> | ||
| 313 | !> @return The value of the derivative | ||
| 314 | 52919796 | pure real(wp) function dbessel_part_dx(xp, yp) | |
| 315 | real(wp), intent(in) :: xp, yp | ||
| 316 | |||
| 317 | dbessel_part_dx = lambda_bessel_mn * bessel_jn_derivative(param_bessel_n, lambda_bessel_mn * xp) * & | ||
| 318 | 52919796 | cos(2 * pi * param_bessel_n * yp) | |
| 319 | 52919796 | end function dbessel_part_dx | |
| 320 | |||
| 321 | !> @brief Compute the derivative of the Bessel part with respect to y | ||
| 322 | !> @param[in] xp The first logical coordinate | ||
| 323 | !> @param[in] yp The second logical coordinate | ||
| 324 | !> | ||
| 325 | !> @return The value of the derivative | ||
| 326 | 52919796 | pure real(wp) function dbessel_part_dy(xp, yp) | |
| 327 | real(wp), intent(in) :: xp, yp | ||
| 328 | |||
| 329 | dbessel_part_dy = -2 * pi * param_bessel_n * bessel_jn(param_bessel_n, lambda_bessel_mn * xp) * & | ||
| 330 | 52919796 | sin(2 * pi * param_bessel_n * yp) | |
| 331 | 52919796 | end function dbessel_part_dy | |
| 332 | |||
| 333 | !> @brief Manufactured solution for the toroidal domain, zero-form | ||
| 334 | !> @param[in] xp The first logical coordinate | ||
| 335 | !> @param[in] yp The second logical coordinate | ||
| 336 | !> @param[in] zp The third logical coordinate | ||
| 337 | !> | ||
| 338 | !> @return The value of the zero-form | ||
| 339 | 22007616 | pure real(wp) function toroidal_zeroform(xp, yp, zp) result(ans) | |
| 340 | real(wp), intent(in) :: xp, yp, zp | ||
| 341 | |||
| 342 | 22007616 | ans = bessel_part(xp, yp) * cos(2 * pi * zp) | |
| 343 | 22007616 | end function toroidal_zeroform | |
| 344 | |||
| 345 | !> @brief Part of the manufactured solution for the toroidal domain, one-form without Bessel part | ||
| 346 | !> @param[in] xp The first logical coordinate | ||
| 347 | !> @param[in] yp The second logical coordinate | ||
| 348 | !> @param[in] zp The third logical coordinate | ||
| 349 | !> | ||
| 350 | !> @return The first component of the one-form without Bessel part | ||
| 351 | 51388572 | pure real(wp) function oneform_x_nobessel(xp, yp, zp) | |
| 352 | real(wp), intent(in) :: xp, yp, zp | ||
| 353 | |||
| 354 | 51388572 | oneform_x_nobessel = cos(2 * pi * yp) * (cos(2 * pi * zp) + sin(2 * pi * zp)) | |
| 355 | 51388572 | end function oneform_x_nobessel | |
| 356 | |||
| 357 | !> @brief Part of the manufactured solution for the toroidal domain, one-form without Bessel part | ||
| 358 | !> @param[in] xp The first logical coordinate | ||
| 359 | !> @param[in] yp The second logical coordinate | ||
| 360 | !> @param[in] zp The third logical coordinate | ||
| 361 | !> | ||
| 362 | !> @return The second component of the one-form without Bessel part | ||
| 363 | 51388572 | pure real(wp) function oneform_y_nobessel(xp, yp, zp) | |
| 364 | real(wp), intent(in) :: xp, yp, zp | ||
| 365 | |||
| 366 | 51388572 | oneform_y_nobessel = -2 * pi * xp * sin(2 * pi * yp) * (cos(2 * pi * zp) + sin(2 * pi * zp)) | |
| 367 | 51388572 | end function oneform_y_nobessel | |
| 368 | |||
| 369 | !> @brief Part of the manufactured solution for the toroidal domain, one-form without Bessel part | ||
| 370 | !> @param[in] xp The first logical coordinate | ||
| 371 | !> @param[in] yp The second logical coordinate | ||
| 372 | !> @param[in] zp The third logical coordinate | ||
| 373 | !> | ||
| 374 | !> @return The third component of the one-form without Bessel part | ||
| 375 | 76151520 | pure real(wp) function oneform_z_nobessel(xp, yp, zp) | |
| 376 | real(wp), intent(in) :: xp, yp, zp | ||
| 377 | |||
| 378 | oneform_z_nobessel = 2 * pi * (coords_major_radius / coords_minor_radius + xp * cos(2 * pi * yp)) *& | ||
| 379 | 76151520 | & (cos(2 * pi * zp) - sin(2 * pi * zp)) | |
| 380 | 76151520 | end function oneform_z_nobessel | |
| 381 | |||
| 382 | !> @brief Manufactured solution for the toroidal domain, one-form | ||
| 383 | !> @param[in] xp The first logical coordinate | ||
| 384 | !> @param[in] yp The second logical coordinate | ||
| 385 | !> @param[in] zp The third logical coordinate | ||
| 386 | !> | ||
| 387 | !> @return The first component of the one-form | ||
| 388 | 24598584 | pure real(wp) function toroidal_oneform_x(xp, yp, zp) result(ans) | |
| 389 | real(wp), intent(in) :: xp, yp, zp | ||
| 390 | |||
| 391 | 24598584 | ans = oneform_x_nobessel(xp, yp, zp) * bessel_part(xp, yp) | |
| 392 | 24598584 | end function toroidal_oneform_x | |
| 393 | |||
| 394 | !> @brief Manufactured solution for the toroidal domain, one-form | ||
| 395 | !> @param[in] xp The first logical coordinate | ||
| 396 | !> @param[in] yp The second logical coordinate | ||
| 397 | !> @param[in] zp The third logical coordinate | ||
| 398 | !> | ||
| 399 | !> @return The second component of the one-form | ||
| 400 | 24598584 | pure real(wp) function toroidal_oneform_y(xp, yp, zp) result(ans) | |
| 401 | real(wp), intent(in) :: xp, yp, zp | ||
| 402 | |||
| 403 | 24598584 | ans = oneform_y_nobessel(xp, yp, zp) * bessel_part(xp, yp) | |
| 404 | 24598584 | end function toroidal_oneform_y | |
| 405 | |||
| 406 | !> @brief Manufactured solution for the toroidal domain, one-form | ||
| 407 | !> @param[in] xp The first logical coordinate | ||
| 408 | !> @param[in] yp The second logical coordinate | ||
| 409 | !> @param[in] zp The third logical coordinate | ||
| 410 | !> | ||
| 411 | !> @return The third component of the one-form | ||
| 412 | 23891904 | pure real(wp) function toroidal_oneform_z(xp, yp, zp) result(ans) | |
| 413 | real(wp), intent(in) :: xp, yp, zp | ||
| 414 | |||
| 415 | 23891904 | ans = oneform_z_nobessel(xp, yp, zp) * bessel_part(xp, yp) | |
| 416 | 23891904 | end function toroidal_oneform_z | |
| 417 | |||
| 418 | !> @brief Manufactured solution for the toroidal domain, two-form | ||
| 419 | !> @param[in] xp The first logical coordinate | ||
| 420 | !> @param[in] yp The second logical coordinate | ||
| 421 | !> @param[in] zp The third logical coordinate | ||
| 422 | !> | ||
| 423 | !> @return The first component of the two-form | ||
| 424 | 26129808 | pure real(wp) function toroidal_twoform_x(xp, yp, zp) result(ans) | |
| 425 | real(wp), intent(in) :: xp, yp, zp | ||
| 426 | |||
| 427 | real(wp) :: dAz_dy, dAy_dz | ||
| 428 | |||
| 429 | 26129808 | dAz_dy = -(2 * pi)**2 * xp * sin(2 * pi * yp) * (cos(2 * pi * zp) - sin(2 * pi * zp)) | |
| 430 | dAy_dz = -(2 * pi)**2 * xp * sin(2 * pi * yp) * (-sin(2 * pi * zp) + cos(2 * pi * zp)) | ||
| 431 | |||
| 432 | ans = bessel_part(xp, yp) * dAz_dy + dbessel_part_dy(xp, yp) * oneform_z_nobessel(xp, yp, zp) & | ||
| 433 | 26129808 | & - bessel_part(xp, yp) * dAy_dz | |
| 434 | 26129808 | end function toroidal_twoform_x | |
| 435 | |||
| 436 | !> @brief Manufactured solution for the toroidal domain, two-form | ||
| 437 | !> @param[in] xp The first logical coordinate | ||
| 438 | !> @param[in] yp The second logical coordinate | ||
| 439 | !> @param[in] zp The third logical coordinate | ||
| 440 | !> | ||
| 441 | !> @return The second component of the two-form | ||
| 442 | 26129808 | pure real(wp) function toroidal_twoform_y(xp, yp, zp) result(ans) | |
| 443 | real(wp), intent(in) :: xp, yp, zp | ||
| 444 | |||
| 445 | real(wp) :: dAx_dz, dAz_dx | ||
| 446 | |||
| 447 | 26129808 | dAx_dz = 2 * pi * cos(2 * pi * yp) * (-sin(2 * pi * zp) + cos(2 * pi * zp)) | |
| 448 | dAz_dx = 2 * pi * cos(2 * pi * yp) * (cos(2 * pi * zp) - sin(2 * pi * zp)) | ||
| 449 | |||
| 450 | ans = bessel_part(xp, yp) * dAx_dz & | ||
| 451 | 26129808 | & - (bessel_part(xp, yp) * dAz_dx + dbessel_part_dx(xp, yp) * oneform_z_nobessel(xp, yp, zp)) | |
| 452 | 26129808 | end function toroidal_twoform_y | |
| 453 | |||
| 454 | !> @brief Manufactured solution for the toroidal domain, two-form | ||
| 455 | !> @param[in] xp The first logical coordinate | ||
| 456 | !> @param[in] yp The second logical coordinate | ||
| 457 | !> @param[in] zp The third logical coordinate | ||
| 458 | !> | ||
| 459 | !> @return The third component of the two-form | ||
| 460 | 26789988 | pure real(wp) function toroidal_twoform_z(xp, yp, zp) result(ans) | |
| 461 | real(wp), intent(in) :: xp, yp, zp | ||
| 462 | |||
| 463 | real(wp) :: dAy_dx, dAx_dy | ||
| 464 | |||
| 465 | 26789988 | dAy_dx = -2 * pi * sin(2 * pi * yp) * (cos(2 * pi * zp) + sin(2 * pi * zp)) | |
| 466 | dAx_dy = -2 * pi * sin(2 * pi * yp) * (cos(2 * pi * zp) + sin(2 * pi * zp)) | ||
| 467 | |||
| 468 | ans = bessel_part(xp, yp) * dAy_dx + dbessel_part_dx(xp, yp) * oneform_y_nobessel(xp, yp, zp) & | ||
| 469 | 26789988 | & - (bessel_part(xp, yp) * dAx_dy + dbessel_part_dy(xp, yp) * oneform_x_nobessel(xp, yp, zp)) | |
| 470 | 26789988 | end function toroidal_twoform_z | |
| 471 | |||
| 472 | !> @brief Manufactured solution for the toroidal domain, three-form | ||
| 473 | !> @param[in] xp The first logical coordinate | ||
| 474 | !> @param[in] yp The second logical coordinate | ||
| 475 | !> @param[in] zp The third logical coordinate | ||
| 476 | !> | ||
| 477 | !> @return The value of the three-form | ||
| 478 | 12608064 | pure real(wp) function toroidal_threeform(xp, yp, zp) result(ans) | |
| 479 | real(wp), intent(in) :: xp, yp, zp | ||
| 480 | |||
| 481 | 12608064 | ans = bessel_part(xp, yp) * cos(2 * pi * zp) * xp * (coords_major_radius / coords_minor_radius + xp * cos(2 * pi * yp)) | |
| 482 | 12608064 | end function toroidal_threeform | |
| 483 | |||
| 484 | end module | ||
| 485 |