Gas model

Project.toml

@@ -10,8 +10,11 @@ GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f"
+PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"

config/simple_gas_valve.yml

@@ -0,0 +1,45 @@
+# Configuration for simple gas valve model
+# This is a sample file. It is not tuned for SPARC.
+# p1 and p2 are taken from DIII-D GASB calibrations where available, and made up
+# when not available. Delay and tau are made up but based (loosely) on experience
+# with the DIII-D gas valve system.
+H2:
+  p1: 6.5268e21  # electrons / s; valve calibration constant
+  p2: 0.61  # V^-1; valve calibration constant
+  delay: 0.0024  # s; time before flow into the VV starts responding to a command
+  tau: 0.071  # s; timescale for flow to catch up to new steady state flow
+D2:
+  p1: 3.325e21  # electrons / s; valve calibration constant
+  p2: 0.78  # V^-1; valve calibration constant
+  delay: 0.0025  # s; time before flow into the VV starts responding to a command
+  tau: 0.072  # s; timescale for flow to catch up to new steady state flow
+N2:
+  p1: 2.628e20
+  p2: 3.41
+  delay: 0.0026
+  tau: 0.073
+CD4:
+  p1: 6.311e20
+  p2: 0.91
+  delay: 0.0025
+  tau: 0.074
+Ne:
+  p1: 1.95e21
+  p2: 0.66
+  delay: 0.0028
+  tau: 0.098
+Ar:
+  p1: 3.992e20
+  p2: 1.0
+  delay: 0.003
+  tau: 0.081
+Kr:
+  p1: 8.75e20
+  p2: 0.825
+  delay: 0.0072
+  tau: 0.121
+Xe:
+  p1: 1.29e20
+  p2: 0.067
+  delay: 0.0103
+  tau: 0.137

src/SD4SOLPS.jl

@@ -11,6 +11,7 @@ export geqdsk_to_imas
 
 include("$(@__DIR__)/supersize_profile.jl")
 include("$(@__DIR__)/repair_eq.jl")
+include("$(@__DIR__)/actuator_model.jl")
 
 greet() = print("Hello World!")
 

src/actuator_model.jl

@@ -0,0 +1,222 @@
+# Actuator models to translate commands (probably in V) into gas flows
+
+using PhysicalConstants.CODATA2018
+using Unitful: Unitful
+using Interpolations: Interpolations
+using YAML: YAML
+
+torrl_per_pam3 = Float64(
+    Unitful.upreferred((1 * Unitful.Torr * Unitful.L) / (Unitful.Pa * Unitful.m^3)),
+)
+electrons_per_molecule = Dict(
+    "H" => 1 * 2,  # Assumed to mean H2. How would you even puff a bunch of H1.
+    "H2" => 1 * 2,
+    "D" => 1 * 2,
+    "D2" => 1 * 2,
+    "T" => 1 * 2,
+    "T2" => 1 * 2,
+    "CH" => 6 + 1,
+    "CD" => 6 + 1,
+    "CH4" => 6 + 4,
+    "CD4" => 6 + 4,
+    "N" => 7 * 2,
+    "N2" => 7 * 2,
+    "Ne" => 10,
+    "Ar" => 18,
+    "Kr" => 36,
+    "Xe" => 54,
+)
+
+"""
+gas_unit_converter()
+
+Converts gas flows between different units. Uses ideal gas law to convert between
+Pressure * volume type flows / quantities and count / current types of units.
+There is a version that accepts floats in and outputs floats, and another that
+deals in Unitful quantities.
+"""
+function gas_unit_converter(
+    value_in::Float64, units_in::String, units_out::String; species::String="H",
+    temperature=293.15,
+)
+    if units_in == units_out
+        return value_in
+    end
+    torrl_to_pam3 = torrl_per_pam3
+
+    pam3_to_molecules = 1.0 / (temperature * BoltzmannConstant.val)
+    torrl_to_molecules = torrl_to_pam3 * pam3_to_molecules
+
+    factor_to_get_molecules_s = Dict(
+        "torr L s^-1" => torrl_to_molecules,
+        "molecules s^-1" => 1.0,
+        "Pa m^3 s^-1" => pam3_to_molecules,
+        "el s^-1" => 1.0 / electrons_per_molecule[species],
+        "A" => ElementaryCharge / electrons_per_molecule[species],
+    )
+    factor_to_get_molecules = Dict(
+        "torr L" => torrl_to_molecules,
+        "molecules" => 1.0,
+        "Pa m^3" => pam3_to_molecules,
+        "el" => 1.0 / electrons_per_molecule[species],
+        "C" => ElementaryCharge / electrons_per_molecule[species],
+    )
+    if haskey(factor_to_get_molecules_s, units_in) &
+       haskey(factor_to_get_molecules_s, units_out)
+        conversion_factor =
+            factor_to_get_molecules_s[units_in] / factor_to_get_molecules_s[units_out]
+    elseif haskey(factor_to_get_molecules, units_in) &
+           haskey(factor_to_get_molecules, units_out)
+        conversion_factor =
+            factor_to_get_molecules[units_in] / factor_to_get_molecules[units_out]
+    else
+        throw(ArgumentError("Unrecognized units: " * units_in * " or " * units_out))
+    end
+    return value_in .* conversion_factor
+end
+
+"""
+gas_unit_converter()
+
+Converts gas flows between different units. Uses ideal gas law to convert between
+Pressure * volume type flows / quantities and count / current types of units.
+This is the Unitful version.
+Output will be unitful, but the units are not simplified automatically. You can
+perform operations such as
+(output |> Unitful.upreferred).val
+Unitful.uconvert(Unitful.whatever, output).val
+to handle simplification or conversion of units.
+
+Although this function pretends torr L s^-1 and Pa m^3 s^-1 are different, use of
+Unitful should cause them to behave the same way as long as you simplify or convert
+units at the end. This means that you can use other pressure*volume type gas units
+and call them torr L s^-1 and the script will deal with them up to having messy
+units in the output.
+"""
+function gas_unit_converter(
+    value_in::Unitful.Quantity, units_in::String, units_out::String;
+    species::String="H", temperature=293.15 * Unitful.K,
+)
+    if units_in == units_out
+        return value_in
+    end
+
+    # The conversion between Torr*L and Pa*m^-3 is the only mundane unit conversion in
+    # here. Unitful or any other unit conversion tool could handle this by itself
+    # (e.g. Unitful.uconvert), but you wouldn't need this fancy function for that.
+    # The rest of the conversions, the hard ones (that require assuming a temperature
+    # or knowing molecule structure and atomic numbers) are all implemented with
+    # multiplicative factors, so Unitful's normal conversion handling is bypassed
+    # and replaced with a multiplicative factor.
+    torrl_to_pam3 = torrl_per_pam3 * Unitful.Pa * Unitful.m^3 / Unitful.Torr / Unitful.L
+    # torrl_to_pam3 should simplify to 1, but Unitful doesn't do this simpliciation
+    # without prompting. So the Torr L or the Pa m^-3 will cancel when this factor
+    # is multiplied by some gas pressure*volume units.
+
+    pam3_to_molecules =
+        Unitful.J / (temperature * BoltzmannConstant) / (Unitful.Pa * Unitful.m^3)
+    torrl_to_molecules = torrl_to_pam3 * pam3_to_molecules
+
+    factor_to_get_molecules_s = Dict(
+        "torr L s^-1" => torrl_to_molecules,
+        "molecules s^-1" => 1.0,
+        "Pa m^3 s^-1" => pam3_to_molecules,
+        "el s^-1" => 1.0 / electrons_per_molecule[species],
+        "A" => ElementaryCharge / electrons_per_molecule[species],
+    )
+    factor_to_get_molecules = Dict(
+        "torr L" => torrl_to_molecules,
+        "molecules" => 1.0,
+        "Pa m^3" => pam3_to_molecules,
+        "el" => 1.0 / electrons_per_molecule[species],
+        "C" => ElementaryCharge / electrons_per_molecule[species],
+    )
+    if haskey(factor_to_get_molecules_s, units_in) &
+       haskey(factor_to_get_molecules_s, units_out)
+        conversion_factor =
+            factor_to_get_molecules_s[units_in] / factor_to_get_molecules_s[units_out]
+    elseif haskey(factor_to_get_molecules, units_in) &
+           haskey(factor_to_get_molecules, units_out)
+        conversion_factor =
+            factor_to_get_molecules[units_in] / factor_to_get_molecules[units_out]
+    else
+        throw(ArgumentError("Unrecognized units: " * units_in * " or " * units_out))
+    end
+    return value_in .* conversion_factor
+end
+
+function select_default_config(model::String)
+    froot = model
+    if model == "instant"
+        froot = "simple"
+    end
+    filename = froot * "_gas_valve.yml"
+    return filename
+end
+
+"""
+    model_gas_valve()
+
+The main function for handling a gas valve model. Has logic for selecting models
+and configurations.
+"""
+function model_gas_valve(
+    model::String; configuration_file::String="auto", species::String="D2",
+)
+    # Select configuration file
+    if configuration_file == "auto"
+        configuration_file = select_default_config(model)
+    end
+    default_config_path = "$(@__DIR__)/../config/"
+    if !occursin("/", configuration_file)
+        configuration_file = default_config_path * configuration_file
+    end
+    if !isfile(configuration_file)
+        throw(ArgumentError("Configuration file not found: " * configuration_file))
+    end
+    config = YAML.load_file(configuration_file)
+    if (species in ["H", "H2", "D", "T", "T2"]) & !(species in keys(config))
+        config = config["D2"]  # They should all be the same
+    elseif (species in ["CH"]) & !(species in keys(config))
+        config = config["CD"]
+    elseif (species in ["CH4", "CD4"]) & !(species in keys(config))
+        config = config["CD4"]
+    else
+        config = config[species]
+    end
+
+    if model == "simple"
+        function simple_gas_model(t, command)
+            flow0 = instant_gas_model(command, config)
+            t_ext = copy(t)
+            prepend!(t_ext, minimum(t) - config["delay"])
+            flow0_ext = copy(flow0)
+            prepend!(flow0_ext, flow0[1])
+            interp = Interpolations.LinearInterpolation(t_ext, flow0_ext)
+            delayed_flow = interp.(t .- config["delay"])
+            return lowpass_filter(t, delayed_flow, config["tau"])
+        end
+        return simple_gas_model
+    elseif model == "instant"
+        instant_gas_model_(t, command) = instant_gas_model(command, config)
+        return instant_gas_model_
+    else
+        throw(ArgumentError("Unrecognized model: " * model))
+    end
+end
+
+function instant_gas_model(command, config)
+    return config["p1"] .* (sqrt.(((command * config["p2"]) .^ 2.0 .+ 1) .- 1))
+end
+
+function lowpass_filter_(raw, previous_smooth, dt, tau)
+    return previous_smooth + dt / (dt + tau) * (raw - previous_smooth)
+end
+
+function lowpass_filter(t, x, tau)
+    xs = zeros(length(t))
+    for i ∈ 2:length(t)
+        xs[i] = lowpass_filter_(x[i], xs[i-1], t[i] - t[i-1], tau)
+    end
+    return xs
+end

test/runtests.jl

@@ -4,6 +4,7 @@ using OMAS: OMAS
 using EFIT: EFIT
 using Plots
 using Test
+using Unitful: Unitful
 
 """
     make_test_profile()
@@ -70,6 +71,51 @@ function define_default_sample_set()
     return b2fgmtry, b2time, b2mn, gridspec, eqdsk
 end
 
+@testset "lightweight_utilities" begin
+    # Gas unit converter
+    flow_tls = 40.63 * Unitful.Torr * Unitful.L / Unitful.s
+    flow_pam3 = SD4SOLPS.gas_unit_converter(flow_tls, "torr L s^-1", "Pa m^3 s^-1")
+    flow_pam3_no_unitful =
+        SD4SOLPS.gas_unit_converter(flow_tls.val, "torr L s^-1", "Pa m^3 s^-1")
+    @test flow_pam3.val > 0.0
+    @test flow_pam3_no_unitful ==
+          Unitful.uconvert(Unitful.Pa * Unitful.m^3 / Unitful.s, flow_pam3).val
+
+    flow_molecules1 = SD4SOLPS.gas_unit_converter(
+        flow_tls,
+        "torr L s^-1",
+        "molecules s^-1";
+        temperature=293.15 * Unitful.K,
+    )
+    flow_molecules2 = SD4SOLPS.gas_unit_converter(
+        flow_tls,
+        "torr L s^-1",
+        "molecules s^-1";
+        temperature=300.0 * Unitful.K,
+    )
+    flow_molecules3 = SD4SOLPS.gas_unit_converter(
+        flow_tls.val,
+        "torr L s^-1",
+        "molecules s^-1";
+        temperature=300.0,
+    )
+    @test flow_molecules1 > flow_molecules2
+    @test (Unitful.upreferred(flow_molecules2)).val == flow_molecules3
+end
+
+@testset "actuator" begin
+    t = collect(LinRange(0, 2.0, 1001))
+    cmd = (t .> 1.0) * 1.55 + (t .> 1.5) * 0.93 + (t .> 1.8) * 0.25
+
+    instant_flow_function = SD4SOLPS.model_gas_valve("instant")
+    instant_flow = instant_flow_function(t, cmd)
+    @test length(instant_flow) == length(cmd)
+
+    simple_flow_function = SD4SOLPS.model_gas_valve("simple")
+    simple_flow = simple_flow_function(t, cmd)
+    @test length(simple_flow) == length(cmd)
+end
+
 @testset "core_profile_extension" begin
     # Just the basic profile extrapolator ------------------
     edge_rho = Array(LinRange(0.88, 1.0, 18))
@@ -140,7 +186,7 @@ end
     psin_out = collect(LinRange(1.0, 1.25, n_outer_prof + 1))[2:end]
 end
 
-@testset "utilities" begin
+@testset "heavy_utilities" begin
     # Test for finding files in allowed folders
     sample_path = splitdir(pathof(SOLPS2IMAS))[1] * "/../samples/"
     file_list = SD4SOLPS.find_files_in_allowed_folders(