diff --git a/src/absorption.f90 b/src/absorption.f90 new file mode 100644 index 0000000..053104f --- /dev/null +++ b/src/absorption.f90 @@ -0,0 +1,400 @@ +! This modules implements a simplified formula for the plasma absorption +! coefficient close to the harmonics valid for low temperature (< 10keV) +! and low harmonic numbers (<6). +! +! Reference: 1983 Nuclear Fusion 23 1153 (https://doi.org/bm6zcr) +! +module absorption + + use const_and_precisions, only : wp_ + + implicit none + + private + public :: alpha_fast + +contains + + pure function line_shape(Y, Npl2, theta) result(phi) + ! The line shape function of the first harmonic, + ! with exception of the X mode in the tenuous plasma regime: + ! + ! φ(ζ) = Im(-1/Z(ζ)) + ! + ! From the same considerations as in `line_shape_harmonics`, we have: + ! + ! φ(ζ) = 1/√π exp(ζ²)/[1 + erfi(ζ)²] + ! + ! Note: Since we already have an implementation of Z(ζ) we + ! use that and write Φ(ζ) = Im(Z)/[Re(Z)² + Im(Z)²]. + ! + use const_and_precisions, only: zero + use dispersion, only: zetac + + implicit none + + ! function arguments + + ! CMA Y variable: Y=ω_c/ω + real(wp_), intent(in) :: Y + ! squared parallel refractive index: N∥²=N²cos²θ + real(wp_), intent(in) :: Npl2 + ! adimensional temperature: Θ = kT/mc² + real(wp_), intent(in) :: theta + ! line shape Φ(ζ) + real(wp_) :: phi + + ! local variables + real(wp_) :: a, b, z + integer :: err + + z = (1 - Y)/sqrt(2 * theta * Npl2) ! ζ + + call zetac(z, zero, a, b, err) ! a=Re(Z(ζ)), b=Im(Z(ζ)) + phi = b/(a**2 + b**2) ! Im(-1/Z) + end function line_shape + + + pure function line_shape_harmonics(Y, Npl2, theta, n) result(phi) + ! The line shape function of all modes at higher harmonics (n>1) + ! and the X mode at the first harmonic: + ! + ! φ(ζ) = Im(Z(ζ)) + ! + ! where Z(ζ) = 1/√π ∫dz exp(-z²)/(z-ζ); + ! ζ = (ω - nω_c)/(√2 k∥ v_t) + ! = (1 - nY)/(√(2Θ)⋅N∥). + ! + ! Since the plasma is weakly dissipative, the wavenumber has a small + ! positive imaginary part, so Im(ζ)<0. By the Sokhotski–Plemelj + ! theorem we have that: + ! + ! lim_Im(ζ)→0 Z(ζ) = 1/√π P ∫dz exp(-z²)/(z-ζ) + i√π exp(-ζ²) + ! = -√π exp(-ζ²) erfi(ζ) + i√π exp(-ζ²) + ! + ! where erfi(ζ) = -i⋅erf(iζ). So, Φ(ζ) = √π exp(-ζ²). + ! + use const_and_precisions, only : sqrt_pi + + implicit none + + ! function arguments + + ! CMA Y variable: Y=ω_c/ω + real(wp_), intent(in) :: Y + ! squared parallel refractive index: N∥²=N²cos²θ + real(wp_), intent(in) :: Npl2 + ! adimensional temperature: Θ = kT/mc² + real(wp_), intent(in) :: theta + ! harmonic number + integer, intent(in) :: n + ! line shape Φ(ζ) + real(wp_) :: phi + + ! local variables + real(wp_) :: z2 + + z2 = (1 - n*Y)**2/(2 * theta * Npl2) + phi = sqrt_pi * exp(-z2) + end function line_shape_harmonics + + + pure function warm_index(d, Npl2, Npr2, sox) result(M2) + ! The real part of the warm plasma refractive index + ! in the finite density regime near the first harmonic: + ! + ! M² = (1 - δ + ½ N⊥²/N² ∓ ½ √Δ) ⋅ N²/N⊥² + ! + ! where Δ = N⊥⁴/N⁴ + 4n²(1-δ)²N∥²/N²; + ! δ = X/Y². + ! + ! Reference: 1983 Nuclear Fusion 23 1153 - table VII + ! + implicit none + + ! function arguments + + ! parameter δ=(ω_pe/ω_ce)² + real(wp_), intent(in) :: d + ! squared refractive index: N∥², N⊥² (cold dispersion) + real(wp_), intent(in) :: Npl2, Npr2 + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! warm refactive indexL M² + real(wp_) :: M2 + + ! local variables + real(wp_) :: N2, delta + + ! total refractive index N² + N2 = Npl2 + Npr2 + + delta = (Npr2/N2)**2 + 4*(1-d)**2 * (Npl2/N2) + M2 = (1 - d + Npr2/N2/2 + sox*sqrt(delta)/2) * N2/Npr2 + end function warm_index + + + pure function warm_index_harmonics(d, Npl2, Npr2, sox, n) result(M2) + ! The real part of the warm plasma refractive index + ! in the finite density regime near the harmonics (n>2) + ! + ! M² = 1 - δ⋅f + ! + ! where f = (1 - δ) / [(1 - δ) - (N⊥²/N² ∓ ρ) / (2n²)]; + ! ρ² = N⊥⁴/N⁴ + 4n²⋅(1 - δ)²⋅ N∥²/N²; + ! δ = X/(nY)². + ! + ! Reference: 1983 Nuclear Fusion 23 1153 - table VIII + ! + implicit none + + ! function arguments + + ! parameter δ=(ω_p/nω_c)² + real(wp_), intent(in) :: d + ! squared refractive index: N∥², N⊥² (cold dispersion) + real(wp_), intent(in) :: Npl2, Npr2 + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! harmonic number + integer, intent(in) :: n + ! warm refactive index + real(wp_) :: M2 + + ! local variables + real(wp_) :: N2, f, r2 + + ! total refractive index N² + N2 = Npl2 + Npr2 + + r2 = (Npr2/N2)**2 + 4*n**2 * (1 - d)**2 * Npl2/N2 + f = (1 - d) / ((1 - d) - (Npr2/N2 + sox * sqrt(r2)) / (2*n**2)) + M2 = 1 - d*f + end function warm_index_harmonics + + + pure function polarisation(X, Y, N2, Npl2, Npr2, sox) result(R) + ! The polarisation function R(θ,O/X) in the finite density regime + ! for the first harmonic: + ! + ! R(N∥, N⊥, O/X) = [(1 - δ)⋅(1 - δ/2 - M∥²) - M⊥²]² + ! ⋅ 1/√(a² + b²) ⋅ M∥ ⋅ 1/d + ! + ! where δ = X/Y² + ! a² = [(1 - δ - M⊥²)² + (1 - d)⋅M∥²]² ⋅ M⊥² + ! b² = (2 - 2δ - M⊥²)² ⋅ (1 - d - M⊥²)² ⋅ M∥² + ! M is the real part of the warm refractive index + ! + ! Reference: 1983 Nuclear Fusion 23 1153 - formula (3.1.68) + ! + implicit none + + ! function arguments + + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: X, Y + ! squared refractive index: N², N∥², N⊥² (cold dispersion) + real(wp_), intent(in) :: N2, Npl2, Npr2 + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! polarisation function + real(wp_) :: R + + ! local variables + real(wp_) :: a2, b2, d + real(wp_) :: M2, Mpl2, Mpr2 + + d = X/Y**2 ! δ=(ω_pe/ω_ce)² + + M2 = warm_index(d, Npl2, Npr2, sox) + Mpl2 = M2 * Npl2/N2 ! M∥² = M²cos²θ + Mpr2 = M2 * Npr2/N2 ! M⊥² = M²sin²θ + + a2 = ((1 - d - Mpr2)**2 + (1 - d)*Mpl2)**2 * Mpr2 + b2 = (2 - 2*d - Mpr2)**2 * (1 - d - Mpr2)**2 * Mpl2 + R = ((1 - d)*(1 - d/2 - Mpl2) - Mpr2)**2 / sqrt(a2 + b2) & + * sqrt(Mpl2) / d + end function polarisation + + + pure function polarisation_harmonics(X, Y, N2, Npl2, Npr2, sox, n) result(mu) + ! The polarisation function μ(n,θ,O/X) times 2(1+cos²θ) + ! in the finite-density plasma regime for harmonics (n>1): + ! + ! μ⋅2(1+cos²θ) = M^(2n-3) (n-1)² (1 - (1+1/n)⋅f) / √(a²+b²) + ! + ! where f as in `warm_index_harmonics`; + ! a = (N⊥/N) (1 + c/e² M∥²); + ! b = (N∥/N) (1 + c/e); + ! c = (1-δ) n²[1 - (1-1/n²)⋅f]²; + ! e = 1-δ-M⊥² + ! M∥ = M/N N∥ = M cosθ; + ! M⊥ = M/N N⊥ = M sinθ; + ! M is the real part of the warm refractive index. + ! + ! Reference: 1983 Nuclear Fusion 23 1153 - formula (3.1.77) + ! + implicit none + + ! function arguments + + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: X, Y + ! squared refractive index: N², N∥², N⊥² (cold dispersion) + real(wp_), intent(in) :: N2, Npl2, Npr2 + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! harmonic number + integer, intent(in) :: n + ! polarisation function + real(wp_) :: mu + + ! local variables + real(wp_) :: d, f, a2, b2, c, e + real(wp_) :: M2, Mpl2, Mpr2 + + ! δ=(ω_pe/nω_ce)² + d = X/(n*Y)**2 + + M2 = warm_index_harmonics(d, Npl2, Npr2, sox, n) + Mpl2 = M2 * Npl2/N2 ! M∥ = Mcosθ + Mpr2 = M2 * Npr2/N2 ! M⊥ = Msinθ + + f = (1 - M2)/d + c = (1 - d) * n**2 * (1 - (1 - 1/n**2) * f)**2 + e = 1 - d - Mpr2 + + a2 = Npr2/N2 * (1 + c/e**2 * Mpl2)**2 ! a² + b2 = Npl2/N2 * (1 + c/e)**2 ! b² + + mu = M2**(n-3.0_wp_/2) * (n-1)**2 & + * (1 - (1 + 1/n) * f)**2 & + / sqrt(a2 + b2) + end function polarisation_harmonics + + + function alpha_fast(X, Y, theta, k0, Npl, Npr, sox, & + nhmin, nhmax) result(alpha) + ! The absorption coefficient around ω=nω_c in the low + ! temperature (< 10keV), low harmonic number limit (< 6), + ! and oblique propagation: + ! + ! α(n,ω,θ,O/X) = α(n,θ)⋅μ(n,θ,O/X)⋅Φ(n,ω) + ! + ! where α(n,θ) is the nth-harmonic absorption coeff.; + ! μ(n,θ,O/X) is the polarisation part; + ! Φ(n,ω) is the line shape function. + ! + ! Finally, for the X1 mode: + ! + ! Notes: + ! 1. The approximation is decent for + ! δ = (ω_pe/nω_ce)² / (2 √Θ N∥) << 1 (tenuous plasma) + ! >> 1 (finite density); + ! 2. In the implementation the angular dependencies + ! have been replaced by cosθ=N∥/N, sinθ=N⊥/N; + ! 3. The special cases for n=1 have been recast + ! from the original form (tbl. XI, eq. 3.1.67) + ! to match with eq. 3.1.73. + ! + implicit none + + ! function arguments + + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: X, Y + ! adimensional temperature: Θ = kT/mc² + real(wp_), intent(in) :: theta + ! vacuum wavenumber k₀=ω/c + real(wp_), intent(in) :: k0 + ! refractive index: N∥, N⊥ (cold dispersion) + real(wp_), intent(in) :: Npl, Npr + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! min/max harmonic number + integer, intent(in) :: nhmin, nhmax + ! absorption coefficient + real(wp_) :: alpha + + ! local variables + integer :: n ! harmonic number + real(wp_) :: N2, Npr2, Npl2 + real(wp_) :: delta + real(wp_) :: shape, polar, absorp + + Npr2 = Npr**2 ! N⊥² refractive index + Npl2 = Npl**2 ! N∥² + N2 = Npr2 + Npl2 ! N² + + ! tenuous/finite density parameter: + ! δ = (ω_p/ω_c)²/(2√Θ N∥) + delta = X/Y**2 / (2*abs(Npl)*sqrt(theta)) + + ! First harmonic + if (delta < 1) then + ! Tenuous plasma + if (sox < 0) then + ! O mode + ! absorption coeff: k₀⋅X + absorp = k0 * X + + ! polarisation part: N⊥⁴/N² ⋅ (N²+2N∥²)²/(N²+ N∥²)³ + polar = Npr2**2/N2 * (N2 + 2*Npl2)**2/(N2 + Npl2)**3 + + ! line shape part: Φ(ζ) = Im(-1/Z(ζ)) / (√(2/Θ)⋅Y⋅|N∥|) + shape = line_shape(Y, Npl2, theta) & + / (sqrt(2/theta) * Y * abs(Npl)) + else + ! X mode + ! polarisation part: (N² + N∥²)/N + polar = (N2 + Npl2) / sqrt(N2) + + ! line shape part: Φ(ζ) = Im(Z(ζ)) / (2√(2θ)⋅Y⋅|N∥|) + shape = line_shape_harmonics(Y, Npl2, theta, n) & + / (2 * sqrt(2*theta) * Y * abs(Npl)) + end if + else + ! Finite density (δ > 1) + + ! absorption coeff: k₀⋅Y + absorp = k0 * Y + + ! polarisation part: R(N∥, N⊥, O/X) + polar = polarisation(X, Y, N2, Npl2, Npr2, sox) + + ! line shape part: Φ(ζ) = Im(-1/Z(ζ)) ⋅ 2√(2Θ) + shape = line_shape(Y, Npl2, theta) & + * 2 * sqrt(2*theta) + end if + alpha = absorp * polar * shape + + if (nhmax == 1) return + + ! n-independent part of α: k₀⋅X + absorp = k0 * X + + ! Higher harmonics + ! note: the loop starts from n=1 to compute + ! the powers/factorials of n, even when nhmin>1. + do n=1,nhmax + ! absorption part: α(n,θ)/[(1+cos²θ)⋅n^(2n)⋅(πω)] + absorp = absorp / (2*n) ! computes (2n)!! + + if (n >= nhmin) then + ! polarisation part: (1+cos²θ)⋅μ(n,θ) + polar = polarisation_harmonics(X, Y, N2, Npl2, Npr2, sox, n) + + ! line shape part: Φ(n,ω)⋅(πω) + shape = line_shape_harmonics(Y, Npl2, theta, n) & + / (sqrt(2*theta) * n*Y * abs(Npl)) + + ! full absorption coefficient: α(n,ω,θ) + alpha = alpha + absorp*n**(2*n) * polar * shape + end if + + ! computes * (Θ N⊥²/N²)^(n-1) + absorp = absorp * (theta * Npr2/N2) + end do + end function alpha_fast + +end module absorption diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 281da95..ff3cba1 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -52,6 +52,7 @@ contains integer, intent(out) :: error ! local variables + ! taucr: max τ before considering a ray fully absorbed real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) character, dimension(2), parameter :: mode=['O','X'] @@ -1985,7 +1986,7 @@ contains real(wp_), intent(in) :: psinv ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: X, Y - ! densityity [10¹⁹ m⁻³], temperature [keV] + ! density [10¹⁹ m⁻³], temperature [keV] real(wp_), intent(in) :: density, temperature ! vacuum wavenumber k₀=ω/c, resonant B field real(wp_), intent(in) :: k0, Bres @@ -2049,28 +2050,48 @@ contains end if end if - ! Compute α from the solution of the dispersion relation - ! The absoption coefficient is defined as - ! - ! α = 2 Im(k̅)⋅s̅ - ! - ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. - ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion - ! relation, we have that - ! - ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| - ! = [2N̅ - ∂(N²s)/∂N∥ b̅ - ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| - ! - ! Assuming Im(k∥)=0: - ! - ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| - ! - block - real(wp_) :: k_im - k_im = k0 * Npr%im ! imaginary part of k⊥ - alpha = 4 * k_im*Npr_cold / derdnm - end block + if (params%iwarm < 4) then + ! Compute α from the solution of the dispersion relation + accurate: block + ! Compute α from the solution of the dispersion relation + ! The absoption coefficient is defined as + ! + ! α = 2 Im(k̅)⋅s̅ + ! + ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. + ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion + ! relation, we have that + ! + ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| + ! = [2N̅ - ∂(N²s)/∂N∥ b̅ + ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| + ! + ! Assuming Im(k∥)=0: + ! + ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| + ! + real(wp_) :: k_im + k_im = k0 * Npr%im ! imaginary part of k⊥ + alpha = 4 * k_im*Npr_cold / derdnm + end block accurate + else + ! Compute a fast and approximate α + fast: block + use logger, only: log_warning + use absorption, only: alpha_fast + + character(256) :: msg + + ! the absorption coefficient in the tenuous plasma limit + alpha = alpha_fast(X, Y, 1/mu, k0, Npl, Npr_cold, sox, nhmin, nhmax) + + if (temperature > 10 .or. nhmax > 5) then + write (msg, '(a,g0.3,a,g0)') & + 'iwarm=4 is inaccurate: Te=', temperature, ' nhmax=', nhmax + call log_warning(msg, mod='gray_core', proc='alpha_effj') + end if + end block fast + end if if (alpha < 0) then error = raise_error(error, negative_absorption)