New version of previous changes

Project.toml

@@ -1,7 +1,7 @@
 name = "HallThruster"
 uuid = "2311f341-5e6d-4941-9e3e-3ce0ae0d9ed6"
 authors = ["Thomas Marks <[email protected]>"]
-version = "0.13.2"
+version = "0.14.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

docs/src/config.md

@@ -40,8 +40,7 @@ Aside from these arguments, all others have  default values provided. These are
 - `source_electron_energy`: Extra source term for electron energy equation. Defaults to `Returns(0.0)`. See [User-Provided Source Terms](@ref) for more information.
 - `LANDMARK`: Whether we are using the LANDMARK physics model. This affects whether certain terms are included in the equations, such as electron and heavy species momentum transfer due to ionization and the form of the electron thermal conductivity. Also affects whether we use an anode sheath model. Defaults to `false`.
 - `ion_wall_losses`: Whether we model ion losses to the walls. Defaults to `false`.
-- `solve_background_neutrals`: Turns on an additional mass flow rate due to neutral ingestion
-- `background_pressure`: The pressure of the background neutrals, in Pascals
+- `background_pressure`: The pressure of the background neutrals, in Pascals. These background neutrals are injected at the anode to simulate the ingestion of facility neutrals.
 - `background_neutral_temperature`: The temperature of the background neutrals, in K
 - `anode_boundary_condition`: Can be either `:sheath` or `:dirichlet`
 - `anom_smoothing_iters`: How many times to smooth the anomalous transport profile

reactions/elastic_Kr.dat

@@ -1,4 +1,5 @@
 Energy (eV)	Rate coefficient (m3/s)
+0.0 0.0
 1.0	1.7652019589294465e-14
 2.0	6.286806105711669e-14
 3.0	1.260621740782443e-13

reactions/elastic_Xe.dat

@@ -1,4 +1,5 @@
 Energy (eV)	Rate coefficient (m3/s)
+0.0 0.0
 1.0	5.130755817456869e-14
 2.0	1.2892031779428643e-13
 3.0	2.2086872081433493e-13

src/HallThruster.jl

@@ -26,7 +26,6 @@ const LANDMARK_FOLDER = joinpath(PACKAGE_ROOT, "landmark")
 const LANDMARK_RATES_FILE = joinpath(LANDMARK_FOLDER, "landmark_rates.csv")
 
 include("utilities/utility_functions.jl")
-include("utilities/transition_functions.jl")
 
 include("physics/physicalconstants.jl")
 include("physics/gas.jl")
@@ -55,14 +54,15 @@ include("simulation/initialization.jl")
 include("simulation/geometry.jl")
 include("simulation/boundaryconditions.jl")
 include("simulation/potential.jl")
-include("simulation/heavy_species.jl")
+include("simulation/update_heavy_species.jl")
 include("simulation/electronenergy.jl")
 include("simulation/sourceterms.jl")
 include("simulation/plume.jl")
-include("simulation/update_electrons.jl")
 include("simulation/configuration.jl")
+include("simulation/update_electrons.jl")
 include("simulation/solution.jl")
 include("simulation/simulation.jl")
+include("simulation/json.jl")
 include("simulation/restart.jl")
 include("simulation/postprocess.jl")
 include("visualization/plotting.jl")
@@ -99,6 +99,7 @@ function example_simulation(;ncells, duration, dt, nsave)
     HallThruster.current_eff(sol_1)
     HallThruster.divergence_eff(sol_1)
     HallThruster.voltage_eff(sol_1)
+    return sol_1
 end
 
 # # Precompile statements to improve load time

src/collisions/anomalous.jl

@@ -58,7 +58,7 @@ function (model::TwoZoneBohm)(νan, params)
         B = params.cache.B[i]
         ωce = e * B / me
 
-        β = params.config.transition_function(params.z_cell[i], params.L_ch, c1, c2)
+        β = linear_transition(params.z_cell[i], params.L_ch, params.config.transition_length, c1, c2)
         νan[i] = β * ωce
     end
 
@@ -221,7 +221,7 @@ function (model::ShiftedTwoZoneBohm)(νan, params)
         B = params.cache.B[i]
         ωce = HallThruster.e * B / HallThruster.me
 
-        β = params.config.transition_function(params.z_cell[i], zstar, c1, c2)
+        β = linear_transition(params.z_cell[i], zstar, params.config.transition_length, c1, c2)
         νan[i] = β * ωce
     end
 end

src/collisions/collision_frequencies.jl

@@ -10,8 +10,8 @@ Effective frequency of electron scattering caused by collisions with neutrals
     return νen
 end
 
-function freq_electron_neutral(U::AbstractArray, params::NamedTuple, i::Int)
-    nn = params.cache.nn_tot[i]
+function freq_electron_neutral(params::NamedTuple, i::Int)
+    nn = params.cache.nn[i]
     Tev = params.cache.Tev[i]
     return freq_electron_neutral(params.electron_neutral_collisions, nn, Tev)
 end

src/numerics/flux.jl

@@ -195,12 +195,10 @@ function compute_edge_states!(UL, UR, U, params; apply_boundary_conditions = fal
 end
 
 function compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions = false)
-    (;config, index, fluids, Δz_edge, num_neutral_fluids, cache) = params
+    (;config, index, fluids, Δz_edge, cache, ncells) = params
     (;λ_global) = cache
-    (;propellant, scheme, ncharge) = config
-    ncells = size(U, 2)
+    (;scheme, ncharge) = config
 
-    mi = propellant.m
     nedges = ncells-1
 
     # Reconstruct the states at the left and right edges using MUSCL scheme
@@ -210,16 +208,14 @@ function compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions = false
     @inbounds for i in 1:nedges
         # Compute wave speeds for each component of the state vector.
         # The only wave speed for neutrals is the neutral convection velocity
-        for j in 1:num_neutral_fluids
-            neutral_fluid = fluids[j]
-            U_neutrals = (U[index.ρn[j], i],)
-            u = velocity(U_neutrals, neutral_fluid)
-            λ_global[j] = abs(u)
-        end
+        neutral_fluid = fluids[1]
+        U_neutrals = (U[index.ρn, i],)
+        u = velocity(U_neutrals, neutral_fluid)
+        λ_global[1] = abs(u)
 
         # Ion wave speeds
         for Z in 1:ncharge
-            fluid_ind = Z + num_neutral_fluids
+            fluid_ind = Z + 1
             fluid = fluids[fluid_ind]
             γ = fluid.species.element.γ
             UL_ions = (UL[index.ρi[Z], i], UL[index.ρiui[Z], i],)
@@ -246,18 +242,16 @@ function compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions = false
 
     @inbounds for i in 1:nedges
         # Neutral fluxes at edge i
-        for j in 1:params.num_neutral_fluids
-            left_state_n  = (UL[index.ρn[j], i],)
-            right_state_n = (UR[index.ρn[j], i],)
+        left_state_n  = (UL[index.ρn, i],)
+        right_state_n = (UR[index.ρn, i],)
 
-            F[index.ρn[j], i] = scheme.flux_function(left_state_n, right_state_n, fluids[j], λ_global[j])[1]
-        end
+        F[index.ρn, i] = scheme.flux_function(left_state_n, right_state_n, fluids[1], λ_global[1])[1]
 
         # Ion fluxes at edge i
         for Z in 1:ncharge
             left_state_i  = (UL[index.ρi[Z], i], UL[index.ρiui[Z], i],)
             right_state_i = (UR[index.ρi[Z], i], UR[index.ρiui[Z], i],)
-            fluid_ind = Z + num_neutral_fluids
+            fluid_ind = Z + 1
             F_mass, F_momentum = scheme.flux_function(left_state_i, right_state_i, fluids[fluid_ind], λ_global[fluid_ind])
             F[index.ρi[Z],   i] = F_mass
             F[index.ρiui[Z], i] = F_momentum

src/numerics/limiters.jl

@@ -17,21 +17,19 @@ const minmod = SlopeLimiter(r -> min(2 / (1 + r), 2r / (1 + r)))
 const koren = SlopeLimiter(r -> max(0, min(2r, min((1 + 2r) / 3, 2))) * 2 / (r+1))
 const osher(β) = SlopeLimiter(r -> (max(0, min(r, β))) * 2 / (r+1))
 
-function stage_limiter!(U, p)
-    ncells = size(U, 2)
-    (;nϵ) = p.cache
-    min_density = p.config.min_number_density * p.config.propellant.m
+function stage_limiter!(U, params)
+    (;ncells, cache, config, index) = params
+    (;nϵ) = cache
+    min_density = config.min_number_density * config.propellant.m
     @inbounds for i in 1:ncells
-        for j in 1:p.num_neutral_fluids
-            U[p.index.ρn[j], j] = max(U[p.index.ρn[j], j], min_density)
-        end
+        U[index.ρn, i] = max(U[index.ρn, i], min_density)
 
-        for Z in 1:p.config.ncharge
-            density_floor = max(U[p.index.ρi[Z], i], min_density)
-            velocity = U[p.index.ρiui[Z], i] / U[p.index.ρi[Z], i]
-            U[p.index.ρi[Z], i] = density_floor
-            U[p.index.ρiui[Z], i] = density_floor * velocity
+        for Z in 1:config.ncharge
+            density_floor = max(U[index.ρi[Z], i], min_density)
+            velocity = U[index.ρiui[Z], i] / U[index.ρi[Z], i]
+            U[index.ρi[Z], i] = density_floor
+            U[index.ρiui[Z], i] = density_floor * velocity
         end
-        nϵ[i] = max(nϵ[i], p.config.min_number_density * p.config.min_electron_temperature)
+        nϵ[i] = max(nϵ[i], config.min_number_density * config.min_electron_temperature)
     end
 end

src/simulation/boundaryconditions.jl

@@ -26,12 +26,10 @@ function left_boundary_state!(bc_state, U, params)
     end
 
     # Add inlet neutral density
-    bc_state[index.ρn[1]] = mdot_a / channel_area[1] / un
+    bc_state[index.ρn] = mdot_a / channel_area[1] / un
 
     # Add ingested mass flow rate at anode
-    if config.solve_background_neutrals
-        bc_state[index.ρn[1]] += params.background_neutral_density * params.background_neutral_velocity / un
-    end
+    bc_state[index.ρn] += params.background_neutral_density * params.background_neutral_velocity / un  * config.neutral_ingestion_multiplier
 
     @inbounds for Z in 1:params.config.ncharge
         interior_density = U[index.ρi[Z],   begin+1]
@@ -63,12 +61,9 @@ function left_boundary_state!(bc_state, U, params)
 
             # Compute boundary flux
             boundary_flux = boundary_velocity * boundary_density
-
-            #boundary_flux = interior_flux
-            #boundary_density = boundary_flux / boundary_velocity
         end
 
-        bc_state[index.ρn[1]] -= boundary_flux / un
+        bc_state[index.ρn] -= boundary_flux / un
         bc_state[index.ρi[Z]] = boundary_density
         bc_state[index.ρiui[Z]] = boundary_flux
     end
@@ -77,10 +72,8 @@ end
 function right_boundary_state!(bc_state, U, params)
     (;index, fluids) = params
 
-    @inbounds for j in 1:params.num_neutral_fluids
-        # Use Neumann boundary conditions for all neutral fluids
-        bc_state[index.ρn[j]] = U[index.ρn[j], end-1]
-    end
+    # Use Neumann boundary conditions for all neutral fluids
+    bc_state[index.ρn] = U[index.ρn, end-1]
 
     @inbounds for Z in 1:params.config.ncharge
         boundary_density = U[index.ρi[Z], end-1]

src/simulation/configuration.jl

@@ -7,7 +7,7 @@ $(TYPEDFIELDS)
 """
 struct Config{A<:AnomalousTransportModel, TC<:ThermalConductivityModel, W<:WallLossModel, IZ<:IonizationModel,
               EX<:ExcitationModel, EN<:ElectronNeutralModel, HET<:Thruster, S_N, S_IC, S_IM, S_ϕ, S_E,
-              T<:TransitionFunction, IC<:InitialCondition, HS<:HyperbolicScheme}
+              IC<:InitialCondition, HS<:HyperbolicScheme}
     discharge_voltage::Float64
     cathode_potential::Float64
     anode_Te::Float64
@@ -27,7 +27,7 @@ struct Config{A<:AnomalousTransportModel, TC<:ThermalConductivityModel, W<:WallL
     electron_ion_collisions::Bool
     min_number_density::Float64
     min_electron_temperature::Float64
-    transition_function::T
+    transition_length::Float64
     initial_condition::IC
     magnetic_field_scale::Float64
     source_neutrals::S_N
@@ -41,9 +41,9 @@ struct Config{A<:AnomalousTransportModel, TC<:ThermalConductivityModel, W<:WallL
     LANDMARK::Bool
     anode_mass_flow_rate::Float64
     ion_wall_losses::Bool
-    solve_background_neutrals::Bool
     background_pressure::Float64
     background_neutral_temperature::Float64
+    neutral_ingestion_multiplier::Float64
     anode_boundary_condition::Symbol
     anom_smoothing_iters::Int
     solve_plume::Bool
@@ -59,7 +59,7 @@ function Config(;
         cathode_potential                   = 0.0,
         cathode_Te                          = 3.0,
         anode_Te                            = cathode_Te,
-        wall_loss_model::WallLossModel      = WallSheath(BNSiO2, 0.15),
+        wall_loss_model::WallLossModel      = WallSheath(BNSiO2, 1.0),
         neutral_velocity                    = 300.0u"m/s",
         neutral_temperature                 = 300.0u"K",
         implicit_energy::Number             = 1.0,
@@ -74,7 +74,7 @@ function Config(;
         electron_ion_collisions::Bool       = true,
         min_number_density                  = 1e6u"m^-3",
         min_electron_temperature            = min(anode_Te, cathode_Te),
-        transition_function::TransitionFunction = LinearTransition(0.2 * thruster.geometry.channel_length, 0.0),
+        transition_length                   = 0.2 * thruster.geometry.channel_length * u"m",
         initial_condition::IC               = DefaultInitialization(),
         magnetic_field_scale::Float64       = 1.0,
         source_neutrals::S_N                = nothing,
@@ -85,9 +85,9 @@ function Config(;
         scheme::HyperbolicScheme            = HyperbolicScheme(),
         LANDMARK                            = false,
         ion_wall_losses                     = false,
-        solve_background_neutrals           = false,
         background_pressure                 = 0.0u"Torr",
         background_neutral_temperature      = 100.0u"K",
+        neutral_ingestion_multiplier::Float64 = 1.0,
         anode_boundary_condition            = :sheath,
         anom_smoothing_iters                = 0,
         solve_plume                         = false,
@@ -101,9 +101,8 @@ function Config(;
     source_IM = ion_source_terms(ncharge, source_ion_momentum,   "momentum")
 
     # Neutral source terms
-    num_neutral_fluids = 1
     if isnothing(source_neutrals)
-        source_neutrals = fill(Returns(0.0), num_neutral_fluids)
+        source_neutrals = fill(Returns(0.0), 1)
     end
 
     # Convert to Float64 if using Unitful
@@ -126,6 +125,8 @@ function Config(;
     background_neutral_temperature = convert_to_float64(background_neutral_temperature, u"K")
     background_pressure = convert_to_float64(background_pressure, u"Pa")
 
+    transition_length = convert_to_float64(transition_length, u"m")
+
     if anode_boundary_condition ∉ [:sheath, :dirichlet, :neumann]
         throw(ArgumentError("Anode boundary condition must be one of :sheath, :dirichlet, or :neumann. Got: $(anode_boundary_condition)"))
     end
@@ -135,7 +136,7 @@ function Config(;
         neutral_velocity, neutral_temperature, implicit_energy, propellant, ncharge,
         ion_temperature, anom_model, conductivity_model,
         ionization_model, excitation_model, electron_neutral_model, electron_ion_collisions,
-        min_number_density, min_electron_temperature, transition_function,
+        min_number_density, min_electron_temperature, transition_length,
         initial_condition, magnetic_field_scale, source_neutrals,
         source_IC,
         source_IM,
@@ -147,9 +148,9 @@ function Config(;
         LANDMARK,
         anode_mass_flow_rate,
         ion_wall_losses,
-        solve_background_neutrals,
         background_pressure,
         background_neutral_temperature,
+        neutral_ingestion_multiplier,
         anode_boundary_condition,
         anom_smoothing_iters,
         solve_plume,
@@ -192,10 +193,10 @@ end
 function configure_fluids(config)
     propellant = config.propellant
 
-    neutral_fluids = [ContinuityOnly(propellant(0); u = config.neutral_velocity, T = config.neutral_temperature)]
+    neutral_fluid = ContinuityOnly(propellant(0); u = config.neutral_velocity, T = config.neutral_temperature)
     ion_fluids = [IsothermalEuler(propellant(Z); T = config.ion_temperature) for Z in 1:config.ncharge]
 
-    fluids = [neutral_fluids; ion_fluids]
+    fluids = [neutral_fluid; ion_fluids]
 
     species = [f.species for f in fluids]
 
@@ -208,7 +209,7 @@ function configure_fluids(config)
 
     last_fluid_index = fluid_ranges[end][end]
     is_velocity_index = fill(false, last_fluid_index)
-    for i in length(neutral_fluids)+2:2:last_fluid_index
+    for i in 3:2:last_fluid_index
         is_velocity_index[i] = true
     end
 
@@ -219,9 +220,7 @@ function configure_index(fluids, fluid_ranges)
     first_ion_fluid_index = findfirst(x -> x.species.Z > 0, fluids)
 
     keys_neutrals = (:ρn, )
-    values_neutrals = (
-        [f[1] for f in fluid_ranges[1:first_ion_fluid_index-1]],
-    )
+    values_neutrals = (1,)
 
     keys_ions = (:ρi, :ρiui)
     values_ions = (