| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module fusion_source | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64, int64 | ||
| 3 | use finite_check, only: is_finite | ||
| 4 | implicit none | ||
| 5 | private | ||
| 6 | public :: dt_reactivity, build_birth_cdf, draw_birth, alpha_energy_joule | ||
| 7 | |||
| 8 | real(dp), parameter :: pi = acos(-1.0_dp) | ||
| 9 | real(dp), parameter :: alpha_energy_joule = 3.5e6_dp*1.602176634e-19_dp | ||
| 10 | |||
| 11 | contains | ||
| 12 | |||
| 13 | 14473 | function dt_reactivity(temperature_ev) result(rate) | |
| 14 | ! H.-S. Bosch & G.M. Hale, Nuclear Fusion 32 (1992) 611-631, | ||
| 15 | ! doi:10.1088/0029-5515/32/4/I07, Eq. (12), Table VII (D-T). | ||
| 16 | ! Input eV; published fit uses keV and cm^3/s; return m^3/s. | ||
| 17 | real(dp), intent(in) :: temperature_ev | ||
| 18 | real(dp) :: rate, t, theta, xi | ||
| 19 | real(dp), parameter :: bg = 34.3827_dp, mrc2 = 1124656.0_dp | ||
| 20 | real(dp), parameter :: c(7) = [1.17302e-9_dp, 1.51361e-2_dp, & | ||
| 21 | 7.51886e-2_dp, 4.60643e-3_dp, 1.35000e-2_dp, & | ||
| 22 | -1.06750e-4_dp, 1.36600e-5_dp] | ||
| 23 | |||
| 24 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14473 times.
|
14473 | if (.not. is_finite(temperature_ev)) & |
| 25 | ✗ | error stop 'fusion source: nonfinite ion temperature' | |
| 26 |
3/4✓ Branch 0 taken 14473 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 14472 times.
|
14473 | if (temperature_ev < 0.0_dp .or. temperature_ev > 1.0e5_dp) & |
| 27 | 1 | error stop 'fusion source: ion temperature outside [0,100000] eV' | |
| 28 | 14472 | rate = 0.0_dp | |
| 29 | 14472 | t = temperature_ev*1.0e-3_dp | |
| 30 | ! Below this temperature exp(-3*xi) underflows to zero. | ||
| 31 |
2/2✓ Branch 0 taken 14464 times.
✓ Branch 1 taken 8 times.
|
14472 | if (t < 1.0e-6_dp) return |
| 32 | theta = t/(1.0_dp - t*(c(2) + t*(c(4) + t*c(6)))/ & | ||
| 33 | 14464 | (1.0_dp + t*(c(3) + t*(c(5) + t*c(7))))) | |
| 34 | 14464 | xi = (bg**2/(4.0_dp*theta))**(1.0_dp/3.0_dp) | |
| 35 | 14464 | rate = 1.0e-6_dp*c(1)*theta*sqrt(xi/(mrc2*t**3))*exp(-3.0_dp*xi) | |
| 36 | 14464 | end function dt_reactivity | |
| 37 | |||
| 38 | 31 | subroutine build_birth_cdf(slo, shi, grid, cdf, power) | |
| 39 | use simple_profiles, only: get_plasma_params | ||
| 40 | use magfie_sub, only: magfie | ||
| 41 | real(dp), intent(in) :: slo, shi | ||
| 42 | integer, intent(in) :: grid(3) | ||
| 43 | real(dp), allocatable, intent(out) :: cdf(:) | ||
| 44 | real(dp), intent(out) :: power | ||
| 45 | real(dp) :: x(3), te, td, tt, nd, nt, source, jac, bmod | ||
| 46 | real(dp) :: bder(3), hcov(3), hctr(3), hcurl(3), cell_volume | ||
| 47 | integer :: ir, it, iz, k, n | ||
| 48 | |||
| 49 |
3/4✓ Branch 0 taken 62 times.
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 62 times.
|
93 | if (.not. all(is_finite([slo, shi]))) & |
| 50 | ✗ | error stop 'fusion source: nonfinite radial bounds' | |
| 51 |
4/6✓ Branch 0 taken 31 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 31 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 30 times.
|
31 | if (slo < 0.0_dp .or. shi > 1.0_dp .or. shi <= slo) & |
| 52 | 1 | error stop 'fusion source: require 0 <= s_inner < s_outer <= 1' | |
| 53 |
9/10✓ Branch 0 taken 88 times.
✓ Branch 1 taken 29 times.
✓ Branch 2 taken 87 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 90 times.
✓ Branch 5 taken 30 times.
✓ Branch 6 taken 90 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 29 times.
|
237 | if (any(grid < 1) .or. any(grid > 1024)) & |
| 54 | 1 | error stop 'fusion source: grid dimensions must be in [1,1024]' | |
| 55 |
3/4✓ Branch 0 taken 87 times.
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 29 times.
|
116 | if (product(int(grid, int64)) > 8388608_int64) & |
| 56 | ✗ | error stop 'fusion source: grid exceeds 8388608 cells' | |
| 57 |
2/2✓ Branch 0 taken 87 times.
✓ Branch 1 taken 29 times.
|
116 | n = product(grid) |
| 58 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 29 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 29 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 29 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 29 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 29 times.
|
87 | allocate (cdf(0:n)) |
| 59 | 29 | cdf(0) = 0.0_dp | |
| 60 | ! Full torus, not one field period. magfie Jacobian is in cm^3. | ||
| 61 | 29 | cell_volume = (shi - slo)*(2.0_dp*pi)**2/real(n, dp)*1.0e-6_dp | |
| 62 | ! Sampling backends read initialized VMEC/Boozer spline tables. Midpoints | ||
| 63 | ! have s > 0, so Boozer's negative-radius flag is never written here. | ||
| 64 | !$omp parallel do default(none) schedule(static) & | ||
| 65 | !$omp shared(slo, shi, grid, cdf, cell_volume, magfie) & | ||
| 66 | !$omp private(ir, it, iz, k, x, te, td, tt, nd, nt, source, & | ||
| 67 | 29 | !$omp jac, bmod, bder, hcov, hctr, hcurl) | |
| 68 | do ir = 1, grid(1) | ||
| 69 | x(1) = slo + (shi - slo)*(real(ir, dp) - 0.5_dp)/grid(1) | ||
| 70 | call get_plasma_params(x(1), te, td, tt, nd, nt) | ||
| 71 | if (.not. all(is_finite([td, tt, nd, nt]))) & | ||
| 72 | error stop 'fusion source: nonfinite profiles' | ||
| 73 | if (min(td, tt, nd, nt) < 0.0_dp) & | ||
| 74 | error stop 'fusion source: negative profiles' | ||
| 75 | ! Relative Maxwellian temperature: mu*(T_D/m_D + T_T/m_T). | ||
| 76 | source = nd*nt*dt_reactivity((3.0_dp*td + 2.0_dp*tt)/5.0_dp) | ||
| 77 | do it = 1, grid(2) | ||
| 78 | x(2) = 2.0_dp*pi*(real(it, dp) - 0.5_dp)/grid(2) | ||
| 79 | do iz = 1, grid(3) | ||
| 80 | x(3) = 2.0_dp*pi*(real(iz, dp) - 0.5_dp)/grid(3) | ||
| 81 | call magfie(x, bmod, jac, bder, hcov, hctr, hcurl) | ||
| 82 | if (.not. is_finite(jac)) & | ||
| 83 | error stop 'fusion source: nonfinite volume Jacobian' | ||
| 84 | k = ((ir - 1)*grid(2) + it - 1)*grid(3) + iz | ||
| 85 | cdf(k) = source*abs(jac)*cell_volume | ||
| 86 | end do | ||
| 87 | end do | ||
| 88 | end do | ||
| 89 | !$omp end parallel do | ||
| 90 | ! Keep the prefix sum ordered: thread count cannot change cell masses | ||
| 91 | ! or their accumulation order, including for fixed-seed sampling. | ||
| 92 |
2/2✓ Branch 0 taken 9256960 times.
✓ Branch 1 taken 29 times.
|
9256989 | do k = 1, n |
| 93 | 9256989 | cdf(k) = cdf(k - 1) + cdf(k) | |
| 94 | end do | ||
| 95 |
3/4✓ Branch 0 taken 29 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 28 times.
|
29 | if (.not. is_finite(cdf(n)) .or. cdf(n) <= 0.0_dp) & |
| 96 | 1 | error stop 'fusion source: total birth rate must be finite and positive' | |
| 97 | 28 | power = cdf(n)*alpha_energy_joule | |
| 98 |
2/2✓ Branch 0 taken 9052185 times.
✓ Branch 1 taken 28 times.
|
9052213 | cdf = cdf/cdf(n) |
| 99 | 28 | end subroutine build_birth_cdf | |
| 100 | |||
| 101 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4400008 times.
|
4400008 | subroutine draw_birth(slo, shi, grid, cdf, z) |
| 102 | real(dp), intent(in) :: slo, shi, cdf(0:) | ||
| 103 | integer, intent(in) :: grid(3) | ||
| 104 | real(dp), intent(out) :: z(5) | ||
| 105 | real(dp) :: u(5) | ||
| 106 | integer :: lo, hi, mid, ir, it, iz, cell | ||
| 107 | |||
| 108 | 4400008 | call random_number(u) | |
| 109 | 4400008 | lo = 0 | |
| 110 | 4400008 | hi = size(cdf) - 1 | |
| 111 |
2/2✓ Branch 0 taken 81049516 times.
✓ Branch 1 taken 4400008 times.
|
85449524 | do while (hi - lo > 1) |
| 112 | 81049516 | mid = (lo + hi)/2 | |
| 113 |
2/2✓ Branch 0 taken 39633842 times.
✓ Branch 1 taken 41415674 times.
|
85449524 | if (cdf(mid) <= u(1)) then |
| 114 | lo = mid | ||
| 115 | else | ||
| 116 | 45815682 | hi = mid | |
| 117 | end if | ||
| 118 | end do | ||
| 119 | 4400008 | cell = hi - 1 | |
| 120 | 4400008 | iz = mod(cell, grid(3)) | |
| 121 | 4400008 | it = mod(cell/grid(3), grid(2)) | |
| 122 | 4400008 | ir = cell/(grid(2)*grid(3)) | |
| 123 | ! Continuous uniform positions inside the selected midpoint cell. | ||
| 124 | 4400008 | z(1) = slo + (shi - slo)*(ir + u(2))/grid(1) | |
| 125 | 4400008 | z(2) = 2.0_dp*pi*(it + u(3))/grid(2) | |
| 126 | 4400008 | z(3) = 2.0_dp*pi*(iz + u(4))/grid(3) | |
| 127 | 4400008 | z(4) = 1.0_dp | |
| 128 | 4400008 | z(5) = 2.0_dp*u(5) - 1.0_dp | |
| 129 | 31 | end subroutine draw_birth | |
| 130 | |||
| 131 | end module fusion_source | ||
| 132 |