Add derived quantities and calculations

src/SynthDiag.jl

@@ -14,4 +14,6 @@ include("$(@__DIR__)/gas_injection.jl")
 
 include("$(@__DIR__)/magic.jl")
 
+include("$(@__DIR__)/derived.jl")
+
 end # module SynthDiag

src/derived.jl

@@ -0,0 +1,124 @@
+# File for holding calculations for derived quantities. The inputs should come from
+# basic inputs/boundary conditions (like injected power) or from synthetic diagnostic
+# outputs (like line average density).
+
+"""
+References
+[Eldon 2022 PPCF] D. Eldon, et al., Plasma Phys. Control. Fusion 64 (2022) 075002 https://doi.org/10.1088/1361-6587/ac6ff9
+[Eich 2013 JNM]
+[Eich 2013 NF]
+"""
+
+export get_powr_from_dd, calc_conducted_loss_power, calc_loss_power, calc_heat_flux_width, calc_q_cyl
+
+function get_power_from_dd(dd::IMASDD.dd)
+    time = dd.summary.time
+    W = dd.summary.global_quantities.energy_mhd.value
+    P_rad_core = dd.summary.global_quantities.power_radiated_inside_lcfs.value
+    P_OHM = dd.summary.global_quantities.power_ohm.value
+    P_NBI =  dd.summary.heating_current_drive.power_nbi.value
+    P_ECH = dd.summary.heating_current_drive.power_ec.value
+    P_ICH = dd.summary.heating_current_drive.power_ic.value
+    P_LH = dd.summary.heating_current_drive.power_lh.value
+    P_fusion = dd.summary.fusion.power.value.value
+    P_other = dd.summary.heating_current_drive.power_additional
+    return time, W, P_rad_core, P_OHM, P_NBI, P_ECH, P_ICH, P_LH, P_fusion, P_other
+end
+
+function set_default_power_arrays(nt::Int64; P_NBI=nothing, P_ECH=nothing, P_ICH=nothing, P_LH=nothing, P_fusion=nothing, P_other=nothing)
+    if P_NBI == nothing
+        P_NBI = zeros(nt)
+    end
+    if P_ECH == nothing
+        P_ECH = zeros(nt)
+    end
+    if P_ICH == nothing
+        P_ICH = zeros(nt)
+    end
+    if P_LH == nothing
+        P_LH = zeros(nt)
+    end
+    if P_fusion == nothing
+        P_fusion = zeros(nt)
+    end
+    if P_other ==  nothing
+        P_other = zeros(nt)
+    end
+    return P_NBI, P_ECH, P_ICH, P_LH, P_fusion, P_other
+end
+
+function calc_conducted_loss_power(
+    time::Array{Float64},
+    W::Array{Float64},
+    P_OHM::Array{Float64};
+    P_NBI=nothing,
+    P_ECH=nothing,
+    P_ICH=nothing,
+    P_LH=nothing,
+    P_fusion=nothing,
+    P_other=nothing,
+)::Array{Float64}
+    nt = length(time)
+    P_NBI, P_ECH, P_ICH, P_LH, P_fusion, P_other = set_default_power_arrays(nt, P_NBI=P_NBI, P_ECH=P_ECH; P_ICH=P_ICH, P_LH=P_LH, P_fusion=P_fusion, P_other=P_other)
+    dWdt = IMASDD.gradient(time, W)
+    P_cond = [calc_conducted_loss_power(dWdt[i], P_OHM[i], P_NBI=P_NBI[i], P_ECH=P_ECH[i], P_ICH=P_ICH[i], P_LH=P_LH[i], P_fusion=P_fusion[i], P_other=P_other[i]) for i in 1:nt]
+    return P_cond
+end
+
+function calc_conducted_loss_power(dWdt::Float64, P_OHM::Float64; P_NBI::Float64=0.0, P_ECH::Float64=0.0, P_ICH::Float64=0.0, P_LH::Float64=0.0, P_fusion::Float64=0.0, P_other::Float64=0.0)::Float64
+    P_input = P_NBI + P_OHM + P_ECH + P_ICH + P_LH + P_fusion + P_other
+    P_cond = P_input - dWdt
+    return P_cond
+end
+
+function calc_loss_power(time::Array{Float64}, W::Array{Float64}, P_rad_core::Array{Float64}, P_OHM::Array{Float64}; P_NBI=nothing, P_ECH=nothing, P_ICH=nothing, P_LH=nothing, P_fusion=nothing, P_other=nothing)::Array{Float64}
+    dWdt = IMASDD.gradient(time, W)
+    nt = length(time)
+    P_NBI, P_ECH, P_ICH, P_LH, P_fusion, P_other = set_default_power_arrays(nt, P_NBI=P_NBI, P_ECH=P_ECH; P_ICH=P_ICH, P_LH=P_LH, P_fusion=P_fusion, P_other=P_other)
+    P_cond = [calc_conducted_loss_power(dWdt[i], P_OHM[i], P_NBI=P_NBI[i], P_ECH=P_ECH[i], P_ICH=P_ICH[i], P_LH=P_LH[i], P_fusion=P_fusion[i], P_other=P_other[i]) for i in 1:nt]
+    P_SOL = P_cond - P_rad_core
+    return P_SOL
+end
+
+function calc_loss_power(dWdt::Float64, P_rad_core::Float64, P_OHM::Float64; P_NBI::Float64=0.0, P_ECH::Float64=0.0, P_ICH::Float64=0.0, P_LH::Float64=0.0, P_fusion::Float64=0.0, P_other::Float64=0.0)::Float64
+    P_cond = calc_conducted_loss_power(dWdt, P_OHM, P_NBI=P_NBI, P_ECH=P_ECH; P_ICH=P_ICH, P_LH=P_LH, P_fusion=P_fusion, P_other=P_other)
+    P_SOL = P_cond - P_rad_core
+    return P_SOL
+end
+
+function calc_loss_power(dd::IMASDD.dd)::Array{Float64}
+    time, W, P_rad_core, P_OHM, P_NBI, P_ECH, P_ICH, P_LH, P_fusion, P_other = get_power_from_dd(dd)
+    P_SOL = calc_loss_power(time, W, P_rad_core, P_OHM, P_NBI=P_NBI, P_ECH=P_ECH; P_ICH=P_ICH, P_LH=P_LH, P_fusion=P_fusion, P_other=P_other)
+    return P_SOL
+end
+
+function calc_q_cyl(B_ϕ_axis::Float64, Iₚ::Float64, aₘᵢₙₒᵣ::Float64, R_geo::Float64, κ::Float64)::Float64
+    # Copied from Equation 19 of [Eldon 2022 PPCF]
+    ε = aₘᵢₙₒᵣ / R_geo
+    μ₀ = π * 4e-7  # H / m
+    q_cyl = π * aₘᵢₙₒᵣ * ε * B_ϕ_axis / (μ₀ * Iₚ) * (1 + κ^2)
+    return q_cyl  # unitless
+end
+
+function calc_heat_flux_width(dd::IMASDD.dd)::Array{Float64}
+    B_ϕ_axis = dd.equilibrium.time_slice[:].global_quantities.magnetic_axis.b_field_tor
+    P_SOL = calc_loss_power(dd)
+    R_geo = dd.summary.boundary.geometric_axis_r.value
+    aₘᵢₙₒᵣ = dd.summary.boundary.minor_radius.value
+    κ = dd.summary.boundary.elongation.value
+    Iₚ = dd.summary.global_quantities.ip.value
+    λq = calc_heat_flux_width.(B_ϕ_axis, P_SOL, R_geo, aₘᵢₙₒᵣ, κ, Iₚ)
+    return λq
+end
+
+function calc_heat_flux_width(B_ϕ_axis::Float64, P_SOL::Float64, R_geo::Float64, aₘᵢₙₒᵣ::Float64, κ::Float64, Iₚ::Float64)::Float64
+    q_cyl = calc_q_cyl(B_ϕ_axis, Iₚ, aₘᵢₙₒᵣ, R_geo, κ)  # Unitless
+    λq = calc_heat_flux_width(B_ϕ_axis, P_SOL, R_geo, q_cyl)
+    return λq
+end
+
+function calc_heat_flux_width(B_ϕ_axis::Float64, P_SOL::Float64, R_geo::Float64, q_cyl::Float64)::Float64
+    # Copied from Equation 18 of [Eldon 2022 PPCF] which came from [Eich 2013 JNM] and [Eich 2013 NF]
+    λq = 0.86 * abs(B_ϕ_axis)^(-0.8) * abs(q_cyl)^(1.11) * abs(P_SOL)^0.11 * abs(R_geo)^(-0.13)
+    return λq
+end

test/runtests.jl

@@ -1,7 +1,8 @@
 using SynthDiag: IMASDD, add_interferometer!, add_langmuir_probes!, add_gas_injection!,
     compute_gas_injection!, get_gas_injection_response, Noise, OverwriteAttemptError,
-    langmuir_probe_current, magic_nesep
-using IMASDD: json2imas
+    langmuir_probe_current, magic_nesep,
+    calc_loss_power, calc_conducted_loss_power, calc_q_cyl, calc_heat_flux_width
+using IMASDD: json2imas, gradient
 using Test
 using Printf
 using Plots
@@ -27,6 +28,9 @@ function parse_commandline()
         ["--magic"],
         Dict(:help => "Test only magic diagnostics",
             :action => :store_true),
+        ["--derived"],
+        Dict(:help => "Test only derived quantities and calculations",
+            :action => :store_true),
     )
     args = ArgParse.parse_args(localARGS, s)
     if !any(values(args)) # If no flags are set, run all tests
@@ -351,3 +355,70 @@ if args["magic"]
         @test length(nesep) == nt
     end
 end
+
+if args["derived"]
+    @testset "derived_quantities" begin
+        # Test data for power calculations
+        time = [i for i in 0:0.01:10.0]
+        Wflat = 2.5e6  # J
+        tau = 0.325  # s
+        ramp_up_end = 1.15
+        ramp_down_start = 9.0
+        W = (time ./ ramp_up_end) .* Wflat
+        W[time .> ramp_up_end] .= Wflat
+        sel = time .> ramp_down_start
+        W[sel] .= Wflat * (1 .- (time[sel] .- ramp_down_start))
+        sel = (time .> 5.5) .& (time .< 6.8)
+        W[sel]  .+= Wflat .* 0.12 .* sin.(time[sel] .* 2 .* π ./ 0.5)
+        dWdt = gradient(time, W)
+        test_time = 0.9
+        test_dWdt = (dWdt[abs.(time .- test_time) .< 0.011])[1]
+
+        Pscale = Wflat / tau
+        P_rad_core = Pscale * 0.32
+        P_input = Pscale + P_rad_core
+        P_NBI = P_input * 0.5
+        P_ECH = P_input * 0.5
+        P_OHM = P_input - (P_NBI + P_ECH)
+
+        nt = length(time)
+
+        P_SOL = calc_loss_power(test_dWdt, P_rad_core, P_OHM, P_NBI=P_NBI, P_ECH=P_ECH)
+        @test P_SOL > 0
+
+        tau_arr = 0.1 .+ time .* 0.0
+        tau_arr[W .> 1e6] .= tau
+        Pscale_arr = W ./ tau_arr
+        P_rad_core_arr = Pscale_arr .* 0.32
+        P_input_arr = Pscale_arr .+ P_rad_core_arr
+        P_NBI_arr = P_input_arr .* 0.5
+        P_ECH_arr = P_input_arr .* 0.4
+        P_OHM_arr = P_input_arr .- (P_NBI_arr .+ P_ECH_arr)
+
+        P_SOL_arr = calc_loss_power(time, W, P_rad_core_arr, P_OHM_arr, P_NBI = P_NBI_arr, P_ECH = P_ECH_arr)
+        @test length(P_SOL_arr) == nt
+
+        P_cond = calc_conducted_loss_power(test_dWdt, P_OHM, P_NBI = P_NBI, P_ECH = P_ECH)
+        @test P_cond > 0
+
+        P_cond_arr = calc_conducted_loss_power(time, W, P_OHM_arr, P_NBI = P_NBI_arr, P_ECH = P_ECH_arr)
+        @test length(P_cond_arr) == nt
+
+        # Test data for qcyl
+        B_ϕ_axis = -2.07  # T
+        Iₚ = 1.1  # MA
+        aₘᵢₙₒᵣ = 0.56  # m
+        R_geo = 1.675  # m
+        κ = 1.81  # unitless
+        q_cyl = calc_q_cyl(B_ϕ_axis, Iₚ, aₘᵢₙₒᵣ, R_geo, κ)
+        @test q_cyl != 0
+
+        λq1 = calc_heat_flux_width(B_ϕ_axis, P_SOL, R_geo, aₘᵢₙₒᵣ, κ, Iₚ)
+        λq2 = calc_heat_flux_width(B_ϕ_axis, P_SOL, R_geo, q_cyl)
+        @test λq1 == λq2
+
+        # Need test for DD version of calc_heat_flux_width
+
+        # Need test for DD version of calc_loss_power
+    end
+end