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,138 @@
+# Actuator models to translate commands (probably in V) into gas flows
+
+using PhysicalConstants.CODATA2018
+using Unitful: Unitful
+using Interpolations: Interpolations
+using YAML: YAML
+
+"""
+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.
+Output will be unitful; take output.val if you do no want the units attached.
+"""
+function gas_unit_converter(
+    value_in, units_in, units_out; species::String="H", temperature=293.15 * Unitful.K,
+)
+    if units_in == units_out
+        return value_in
+    end
+
+    # Multiply by gas flow to convert Torr L/s to Pa m^3/s
+    torrl_to_pam3 = 0.133322368 * Unitful.Pa * Unitful.m^3 / Unitful.Torr / Unitful.L
+    pam3_to_molecules =
+        Unitful.J / (temperature * BoltzmannConstant) / (Unitful.Pa * Unitful.m^3)
+    torrl_to_molecules = torrl_to_pam3 * pam3_to_molecules
+    atoms_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,
+    )
+    atoms_per_molecule[species]
+
+    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 / atoms_per_molecule[species],
+        "A" => ElementaryCharge / atoms_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 / atoms_per_molecule[species],
+        "C" => ElementaryCharge / atoms_per_molecule[species],
+    )
+    if haskey(factor_to_get_molecules_s, units_in)
+        conversion_factor =
+            factor_to_get_molecules_s[units_in] / factor_to_get_molecules_s[units_out]
+    elseif haskey(factor_to_get_molecules, units_in)
+        conversion_factor =
+            factor_to_get_molecules[units_in] / factor_to_get_molecules[units_out]
+    else
+        throw(ArgumentError("Unrecognized units: " * units_in))
+    end
+    return value_in .* conversion_factor
+end
+
+select_default_config(model::String) = model * "_gas_valve.yml"
+
+"""
+    model_gas_valve()
+
+The main function for handling a gas valve model. Has logic for selecting models
+and configurations.
+"""
+function model_gas_valve(
+    t, command, 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"
+        return simple_gas_model(t, command, config)
+    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
+
+function simple_gas_model(t, command, config)
+    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 flow = lowpass_filter(t, delayed_flow, config["tau"])
+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,33 @@ 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")
+    @test flow_pam3.val > 0.0
+    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,
+    )
+    @test flow_molecules1 > flow_molecules2
+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
+    flow = SD4SOLPS.model_gas_valve(t, cmd, "simple")
+    @test length(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 +168,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(