Reactions refactor

docs/src/collisions.md

@@ -212,26 +212,30 @@ This method returns an integer corresponding to the maximum allowed charge state
 ```julia
 import HallThruster: maximum_charge_state
 
-maximum_charge_state(::MyIonizationModel) = 1
+HallThruster.maximum_charge_state(::MyIonizationModel) = 1
 ```
 
 ### `load_reactions(model::ReactionModel, species::Vector{Species})`
 
-This is the most important method to define. `load_reactions` takes our model along with a vector of `Species` (generated using the `propellant` and `ncharge` fields in the user-provided `Config`) and returns a vector of `IonizationReaction{I}` or `ExcitationReaction{I}`, where `I` is the type of the rate coefficient function. It is important that the `Reactions` all have the same type so that the returned vector will have a concrete type, which will prevent dynamic dispatch during runtime and prevent unnecessary slowdowns. This means that if you have multiple reactions with different rate coefficient functions, you should use something like [`FunctionWrappers.jl`](https://github.com/yuyichao/FunctionWrappers.jl).  Let's implement a simple ionization curve with the form
+`load_reactions` takes our model along with a vector of `Species` (generated using the `propellant` and `ncharge` fields in the user-provided `Config`) and returns a vector of `IonizationReaction` or `ExcitationReaction`. If you are not using a lookup table for this purpose, the `rate_coeffs` field
+of these structs can be left empty. In that case you also need to define the `rate_coeff` function, as by default that looks for a lookup table in the reaction struct.
 
+Let's say our ionization rate coefficient is
 ```math
 k_{iz}(\epsilon) = 4\times 10^{-20} \exp\left(\frac{-7.3}{\frac{2}{3} \epsilon}\right)\sqrt{\frac{8 (\frac{2}{3}\epsilon)}{\pi m_e}}
 ```
 
 We would implement this as so:
 
 ```julia
-using FunctionWrappers
-import HallThruster: e, me, load_reactions
+import HallThruster: e, me, load_reactions, rate_coefficient
 
-kiz(ϵ) = 4e-20 * exp(-7.3 / (2/3 * ϵ)) * sqrt(8 * 2/3 * ϵ / pi / me)
+function rate_coeff(::MyIonizationModel, ::Reaction, energy::Float64)
+    Tev = 2/3 * energy
+    return 4e-20 * exp(-7.3 / Tev) * sqrt(8 * Tev / pi / me)
+end
 
-function load_reactions(model::MyIonizationModel, species)
+function load_reactions(::MyIonizationModel, species)
     rxn = IonizationReaction(
     	energy = -7.3,
         #=
@@ -241,8 +245,8 @@ function load_reactions(model::MyIonizationModel, species)
         =#
         reactant = species[1],
         product = species[2],
-        # Use a function wrapper here, though not necessary with only one reaction
-        rate_coeff = FunctionWrapper{Float64, Tuple{Float64}}(kiz)
+        # Empty lookup table, as we will not be using it
+        rate_coeff = Float64[]
     )
     return rxn
 end

src/collisions/collision_frequencies.jl

@@ -2,18 +2,18 @@
     freq_electron_neutral(model::ElectronNeutralModel, nn, Tev)
 Effective frequency of electron scattering caused by collisions with neutrals
 """
-@inline function freq_electron_neutral(collisions::Vector{ElasticCollision{T}}, nn::Number, Tev::Number) where T
+@inline function freq_electron_neutral(model::ElectronNeutralModel, collisions::Vector{ElasticCollision}, nn::Number, Tev::Number)
     νen = 0.0
     @inbounds for c in collisions
-        νen += c.rate_coeff(3/2 * Tev) * nn
+        νen += rate_coeff(model, c, 3/2 * Tev) * nn
     end
     return νen
 end
 
-function freq_electron_neutral(params::NamedTuple, i::Int)
+function freq_electron_neutral(params, i)
     nn = params.cache.nn[i]
     Tev = params.cache.Tev[i]
-    return freq_electron_neutral(params.electron_neutral_collisions, nn, Tev)
+    return freq_electron_neutral(params.config.electron_neutral_model, params.electron_neutral_collisions, nn, Tev)
 end
 
 """

src/collisions/elastic.jl

@@ -1,48 +1,39 @@
-Base.@kwdef struct ElasticCollision{I} <: Reaction
+Base.@kwdef struct ElasticCollision <: Reaction
     species::Species
-    rate_coeff::I
+    rate_coeffs::Vector{Float64}
 end
 
 abstract type ElectronNeutralModel <: ReactionModel end
 
+# Definitions for NoElectronNeutral
 struct NoElectronNeutral <: ElectronNeutralModel end
-
 supported_species(::NoElectronNeutral) = Gas[]
+load_reactions(::NoElectronNeutral, species) = ElasticCollision[]
+rate_coeff(::NoElectronNeutral, ::Reaction, ::Float64) = 0.0
 
-load_reactions(::NoElectronNeutral, species) = ElasticCollision{Nothing}[]
-
-
+# Definitions for LandmarkElectronNeutral
 struct LandmarkElectronNeutral <: ElectronNeutralModel end
-
 supported_species(::LandmarkElectronNeutral) = [Xenon]
+load_reactions(::LandmarkElectronNeutral, species) = [ElasticCollision(Xenon(0), Float64[])]
+rate_coeff(::LandmarkElectronNeutral, ::Reaction, ::Float64) = 2.5e-13
 
-function load_reactions(::LandmarkElectronNeutral, species)
-    landmark_electron_neutral = Returns(2.5e-13)
-    return [ElasticCollision(Xenon(0), landmark_electron_neutral)]
-end
-
-
-struct GKElectronNeutral <: ElectronNeutralModel end
-
-supported_species(::GKElectronNeutral) = [Xenon]
-
+# Definitions for GKElectronNeutral
 """
-    σ_en(Tev)
-
-Electron neutral collision cross section in m² as a function of electron temperature in eV.
+Electron neutral collision model as defined by
 Eq. 3.6-13, from Fundamentals of Electric Propulsion, Goebel and Katz, 2008.
-
 """
-@inline function σ_en(Tev)
-    return max(0.0, 6.6e-19 * ((Tev / 4 - 0.1) / (1 + (Tev / 4)^1.6)))
-end
+struct GKElectronNeutral <: ElectronNeutralModel end
+supported_species(::GKElectronNeutral) = [Xenon]
+load_reactions(::GKElectronNeutral, species) = [ElasticCollision(Xenon(0), Float64[])]
 
-function load_reactions(::GKElectronNeutral, species)
-    rate_coeff(ϵ) = σ_en(2/3 * ϵ) * electron_sound_speed(2/3 * ϵ)
-    return [ElasticCollision(Xenon(0), rate_coeff)]
-end
+@inline σ_en(Tev) = max(0.0, 6.6e-19 * ((Tev / 4 - 0.1) / (1 + (Tev / 4)^1.6)))
 
+@inline function rate_coeff(::GKElectronNeutral, ::Reaction, energy::Float64)
+    Tev = 2/3 * energy
+    return σ_en(Tev) * electron_sound_speed(Tev)
+end
 
+# Definitions for ElectronNeutralLookup
 Base.@kwdef struct ElectronNeutralLookup <: ElectronNeutralModel
     directories::Vector{String} = String[]
 end
@@ -59,15 +50,14 @@ function load_reactions(model::ElectronNeutralLookup, species)
         reactant = species_sorted[i]
         if reactant.Z > 0
             break
-        else
-            for folder in folders
-                filename = rate_coeff_filename(reactant, product, collision_type, folder)
-                if ispath(filename)
-                    _, rate_coeff = load_rate_coeffs(reactant, product, collision_type,folder)
-                    reaction = HallThruster.ElasticCollision(reactant, rate_coeff)
-                    push!(reactions, reaction)
-                    break
-                end
+        end
+        for folder in folders
+            filename = rate_coeff_filename(reactant, product, collision_type, folder)
+            if ispath(filename)
+                _, rate_coeff = load_rate_coeffs(reactant, product, collision_type,folder)
+                reaction = HallThruster.ElasticCollision(reactant, rate_coeff)
+                push!(reactions, reaction)
+                break
             end
         end
     end

src/collisions/excitation.jl

@@ -1,7 +1,7 @@
-Base.@kwdef struct ExcitationReaction{I} <: Reaction
+Base.@kwdef struct ExcitationReaction <: Reaction
     energy::Float64
     reactant::Species
-    rate_coeff::I
+    rate_coeffs::Vector{Float64}
 end
 
 function Base.show(io::IO, i::ExcitationReaction)
@@ -64,7 +64,7 @@ Supports only singly-charged Xenon.
 """
 struct LandmarkExcitationLookup <: ExcitationModel end
 
-function load_reactions(model::LandmarkExcitationLookup, species)
+function load_reactions(::LandmarkExcitationLookup, species)
     rates = readdlm(LANDMARK_RATES_FILE, ',', skipstart = 1)
     ϵ = rates[:, 1]
     k_iz = rates[:, 2]
@@ -75,6 +75,8 @@ function load_reactions(model::LandmarkExcitationLookup, species)
     # We infer the excitation loss coefficient by subtracting the contribution of the ionization
     # loss coefficient from the inelastic loss term
     k_excitation = @. (k_loss - ionization_energy * k_iz) / excitation_energy
-    rate_coeff = LinearInterpolation(ϵ, k_excitation)
-    return [ExcitationReaction(excitation_energy, Xenon(0), rate_coeff)]
+    itp = LinearInterpolation(ϵ, k_excitation)
+    xs = 0:1.0:255
+    rate_coeffs = itp.(xs)
+    return [ExcitationReaction(excitation_energy, Xenon(0), rate_coeffs)]
 end

src/collisions/ionization.jl

@@ -1,8 +1,8 @@
-Base.@kwdef struct IonizationReaction{I} <: Reaction
+Base.@kwdef struct IonizationReaction <: Reaction
     energy::Float64
     reactant::Species
     product::Species
-    rate_coeff::I
+    rate_coeffs::Vector{Float64}
 end
 
 function Base.show(io::IO, i::IonizationReaction)
@@ -74,6 +74,8 @@ function load_reactions(::LandmarkIonizationLookup, species)
     rates = readdlm(LANDMARK_RATES_FILE, ',', skipstart = 1)
     ϵ = rates[:, 1]
     k = rates[:, 2]
-    rate_coeff = LinearInterpolation(ϵ, k)
-    return [IonizationReaction(12.12, Xenon(0), Xenon(1), rate_coeff)]
+    itp = LinearInterpolation(ϵ, k)
+    xs = 0:1.0:255
+    rate_coeffs = itp.(xs)
+    return [IonizationReaction(12.12, Xenon(0), Xenon(1), rate_coeffs)]
 end

src/collisions/reactions.jl

@@ -1,5 +1,6 @@
 abstract type Reaction end
 
+
 function rate_coeff_filename(reactant, product, reaction_type, folder = REACTION_FOLDER)
     fname = if product === nothing
         joinpath(folder, join([reaction_type, repr(reactant)], "_") * ".dat")
@@ -56,12 +57,27 @@ function load_rate_coeffs(reactant, product, reaction_type, folder = REACTION_FO
     end
     ϵ = rates[:, 1]
     k = rates[:, 2]
-    rate_coeff = LinearInterpolation(ϵ, k, resample_uniform=true)
-    return energy, rate_coeff
+    itp = LinearInterpolation(ϵ, k)
+    xs = 0:1.0:255
+    rate_coeffs = itp.(xs)
+    return energy, rate_coeffs
 end
 
 abstract type ReactionModel end
 
+"""
+By default, rate_coeff looks for a lookup table stored in the reaction struct
+"""
+function rate_coeff(::ReactionModel, rxn::Reaction, energy)
+    ind = Base.trunc(Int64, energy)
+    N = length(rxn.rate_coeffs) - 2
+    ind = ind > N ? N : ind
+    r1 = rxn.rate_coeffs[ind+1]
+    r2 = rxn.rate_coeffs[ind+2]
+    t = energy - ind
+    return (1 - t) * r1 + t * r2
+end
+
 """
     load_reactions(model::ReactionModel, species)::Vector{IonizationReaction}
 Load ionization reactions for the provided species and ionization model
@@ -109,13 +125,11 @@ function _load_reactions(model::ReactionModel, species)
 end
 
 function _indices(symbol, reactions, species_range_dict)
-    indices = [Int[] for _ in reactions]
+    indices = zeros(Int, length(reactions))
     for (i, reaction) in enumerate(reactions)
         species = getfield(reaction, symbol).symbol
-        ranges = species_range_dict[species]
-        for r in ranges
-            push!(indices[i], r[1])
-        end
+        range = species_range_dict[species]
+        indices[i] = range[1]
     end
     return indices
 end

src/simulation/configuration.jl

@@ -111,7 +111,7 @@ function Config(;
 
     anode_Te = convert_to_float64(anode_Te, u"eV")
     cathode_Te = convert_to_float64(cathode_Te, u"eV")
-    
+
     default_neutral_velocity = 150.0 # m/s
     default_neutral_temp = 500.0 # K
     if isnothing(neutral_velocity) && isnothing(neutral_temperature)
@@ -217,10 +217,10 @@ function configure_fluids(config)
     species = [f.species for f in fluids]
 
     fluid_ranges = ranges(fluids)
-    species_range_dict = Dict(Symbol(fluid.species) => UnitRange{Int64}[] for fluid in fluids)
+    species_range_dict = Dict(Symbol(fluid.species) => 0:0 for fluid in fluids)
 
     for (fluid, fluid_range) in zip(fluids, fluid_ranges)
-        push!(species_range_dict[Symbol(fluid.species)], fluid_range)
+        species_range_dict[Symbol(fluid.species)] = fluid_range
     end
 
     last_fluid_index = fluid_ranges[end][end]

src/simulation/electronenergy.jl

@@ -22,10 +22,13 @@ function update_electron_energy!(params, dt)
 
     bϵ[end] = 1.5 * Te_R * ne[end]
 
+    # Compute energy source terms
+    Q = cache.cell_cache_1
+    source_electron_energy!(Q, params)
+
     @inbounds for i in 2:ncells-1
-        Q = source_electron_energy(params, i)
         # User-provided source term
-        Q += config.source_energy(params, i)
+        Q[i] += config.source_energy(params, i)
 
         neL = ne[i-1]
         ne0 = ne[i]
@@ -117,7 +120,7 @@ function update_electron_energy!(params, dt)
         flux = 5/3 * nϵ0 * ue0
 
         # Explicit right-hand-side
-        bϵ[i] = nϵ[i] + dt * (Q - explicit * F_explicit)
+        bϵ[i] = nϵ[i] + dt * (Q[i] - explicit * F_explicit)
         bϵ[i] -= dt * flux * dlnA_dz
     end
 

src/simulation/restart.jl

@@ -34,14 +34,14 @@ function read_restart(path::AbstractString)
     )
 end
 
-function load_restart(grid, fluids, config, path::AbstractString)
+function load_restart(grid, config, path::AbstractString)
     sol = read_restart(path)
-    U, cache = load_restart(grid, fluids, config, sol)
+    U, cache = load_restart(grid, config, sol)
     return U, cache
 end
 
-function load_restart(grid, fluids, config, sol::Solution)
-    U, cache = HallThruster.allocate_arrays(grid, fluids)
+function load_restart(grid, config, sol::Solution)
+    U, cache = HallThruster.allocate_arrays(grid, config)
 
     z_cell = sol.params.z_cell