From 6108a00db1ee714c916316f6b9195381b7fc2f53 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 20:05:03 -0400 Subject: [PATCH 01/14] remove multiple neutral fluids --- HallThrusterTutorial.ipynb | 6 +++--- src/collisions/collision_frequencies.jl | 2 +- src/numerics/flux.jl | 24 +++++++++------------ src/numerics/limiters.jl | 4 +--- src/simulation/boundaryconditions.jl | 12 +++++------ src/simulation/configuration.jl | 13 +++++------ src/simulation/heavy_species.jl | 11 ++++------ src/simulation/initialization.jl | 2 +- src/simulation/postprocess.jl | 7 ++---- src/simulation/simulation.jl | 8 ++----- src/simulation/solution.jl | 4 ++-- src/simulation/sourceterms.jl | 17 +++++---------- src/simulation/update_electrons.jl | 14 +++++------- src/visualization/recipes.jl | 4 ++-- test/order_verification/ovs_ions.jl | 5 ++--- test/unit_tests/test_boundary_conditions.jl | 9 ++++---- test/unit_tests/test_electrons.jl | 2 +- test/unit_tests/test_initialization.jl | 2 +- test/unit_tests/test_misc.jl | 16 +++++++------- test/unit_tests/test_walls.jl | 6 +++--- 20 files changed, 67 insertions(+), 101 deletions(-) diff --git a/HallThrusterTutorial.ipynb b/HallThrusterTutorial.ipynb index d6d52c00..a6806560 100644 --- a/HallThrusterTutorial.ipynb +++ b/HallThrusterTutorial.ipynb @@ -5029,15 +5029,15 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.10.0", + "display_name": "Julia 1.11.0-beta1", "language": "julia", - "name": "julia-1.10" + "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.10.0" + "version": "1.11.0-beta1" } }, "nbformat": 4, diff --git a/src/collisions/collision_frequencies.jl b/src/collisions/collision_frequencies.jl index f4addb11..73c084b5 100644 --- a/src/collisions/collision_frequencies.jl +++ b/src/collisions/collision_frequencies.jl @@ -11,7 +11,7 @@ Effective frequency of electron scattering caused by collisions with neutrals end function freq_electron_neutral(U::AbstractArray, params::NamedTuple, i::Int) - nn = params.cache.nn_tot[i] + nn = params.cache.nn[i] Tev = params.cache.Tev[i] return freq_electron_neutral(params.electron_neutral_collisions, nn, Tev) end diff --git a/src/numerics/flux.jl b/src/numerics/flux.jl index 0158f7f6..860ba879 100644 --- a/src/numerics/flux.jl +++ b/src/numerics/flux.jl @@ -195,7 +195,7 @@ 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) = params (;λ_global) = cache (;propellant, scheme, ncharge) = config ncells = size(U, 2) @@ -210,16 +210,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 +244,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 diff --git a/src/numerics/limiters.jl b/src/numerics/limiters.jl index de860435..2e935b05 100644 --- a/src/numerics/limiters.jl +++ b/src/numerics/limiters.jl @@ -22,9 +22,7 @@ function stage_limiter!(U, p) (;nϵ) = p.cache min_density = p.config.min_number_density * p.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[p.index.ρn, i] = max(U[p.index.ρn, i], min_density) for Z in 1:p.config.ncharge density_floor = max(U[p.index.ρi[Z], i], min_density) diff --git a/src/simulation/boundaryconditions.jl b/src/simulation/boundaryconditions.jl index 2b441a9a..d35bf608 100644 --- a/src/simulation/boundaryconditions.jl +++ b/src/simulation/boundaryconditions.jl @@ -26,11 +26,11 @@ 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 + bc_state[index.ρn] += params.background_neutral_density * params.background_neutral_velocity / un end @inbounds for Z in 1:params.config.ncharge @@ -68,7 +68,7 @@ function left_boundary_state!(bc_state, U, params) #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 +77,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] diff --git a/src/simulation/configuration.jl b/src/simulation/configuration.jl index 178ef831..0ae1f964 100644 --- a/src/simulation/configuration.jl +++ b/src/simulation/configuration.jl @@ -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 @@ -192,10 +191,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 +207,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 +218,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 = ( diff --git a/src/simulation/heavy_species.jl b/src/simulation/heavy_species.jl index e9cc36e9..5725db23 100644 --- a/src/simulation/heavy_species.jl +++ b/src/simulation/heavy_species.jl @@ -23,14 +23,11 @@ function update_heavy_species!(dU, U, params, t; apply_boundary_conditions = tru dlnA_dz = dA_dz[i] / channel_area[i] - # Handle neutrals - for j in 1:params.num_neutral_fluids - # Neutral fluxes - dU[index.ρn[j], i] = (F[index.ρn[j], left] - F[index.ρn[j], right]) / Δz + # Neutral flux + dU[index.ρn, i] = (F[index.ρn, left] - F[index.ρn, right]) / Δz - # User-provided neutral source term - dU[index.ρn[j], i] += source_neutrals[j](U, params, i) - end + # User-provided neutral source term + dU[index.ρn, i] += source_neutrals[1](U, params, i) # Handle ions for Z in 1:ncharge diff --git a/src/simulation/initialization.jl b/src/simulation/initialization.jl index 719bb65f..7ea37476 100644 --- a/src/simulation/initialization.jl +++ b/src/simulation/initialization.jl @@ -58,7 +58,7 @@ function initialize!(U, params, ::DefaultInitialization) # Fill the state vector for (i, z) in enumerate(z_cell) - U[index.ρn[1], i] = neutral_function(z) + U[index.ρn, i] = neutral_function(z) for Z in 1:params.config.ncharge U[index.ρi[Z], i] = ion_density_function(z, Z) diff --git a/src/simulation/postprocess.jl b/src/simulation/postprocess.jl index d8ff96d0..169feebd 100644 --- a/src/simulation/postprocess.jl +++ b/src/simulation/postprocess.jl @@ -158,9 +158,7 @@ function Base.getindex(sol::Solution, field::Symbol, charge::Int = 1) throw(ArgumentError("No ions of charge state $charge in Hall thruster solution. Maximum charge state in provided solution is $(sol.params.config.ncharge).")) end - if field == :nn - return [saved[:nn][charge, :] for saved in sol.savevals] - elseif field == :ni + if field == :ni return [saved[:ni][charge, :] for saved in sol.savevals] elseif field == :ui return [saved[:ui][charge, :] for saved in sol.savevals] @@ -225,8 +223,7 @@ function load_landmark_data(case, suffix; ncells = 100) niui = ui' |> collect, ∇ϕ = -E_itp.(zs), ϕ = ϕ_itp.(zs), - nn_tot = nns, - nn = nns' |> collect, + nn = nns, ) ionization_reactions = HallThruster.load_reactions(LandmarkIonizationLookup(), [Xenon(0), Xenon(1)]); diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index f8d0e94a..39685a1f 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -50,13 +50,11 @@ function allocate_arrays(grid, fluids, anom_model = HallThruster.NoAnom()) wall_losses = zeros(ncells) inelastic_losses = zeros(ncells) - num_neutral_fluids = count(f -> f.species.Z == 0, fluids) ncharge = maximum(f.species.Z for f in fluids) ni = zeros(ncharge, ncells) ui = zeros(ncharge, ncells) niui = zeros(ncharge, ncells) - nn = zeros(num_neutral_fluids, ncells) - nn_tot = zeros(ncells) + nn = zeros(ncells) γ_SEE = zeros(ncells) Id = [0.0] error_integral = [0.0] @@ -95,7 +93,7 @@ function allocate_arrays(grid, fluids, anom_model = HallThruster.NoAnom()) cache = (; Aϵ, bϵ, nϵ, B, νan, νc, μ, ϕ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, νiw, νe, κ, F, UL, UR, Z_eff, λ_global, νiz, νex, K, Id, ji, - ni, ui, Vs, niui, nn, nn_tot, k, u1, γ_SEE, cell_cache_1, + ni, ui, Vs, niui, nn, k, u1, γ_SEE, cell_cache_1, error_integral, Id_smoothed, anom_multiplier, smoothing_time_constant, errors, dcoeffs, ohmic_heating, wall_losses, inelastic_losses, @@ -144,7 +142,6 @@ function run_simulation( dt = convert_to_float64(dt, u"s") fluids, fluid_ranges, species, species_range_dict, is_velocity_index = configure_fluids(config) - num_neutral_fluids = count(f -> f.species.Z == 0, fluids) if (grid.ncells == 0) # backwards compatability for unspecified grid type @@ -230,7 +227,6 @@ function run_simulation( dt = [dt], CFL, adaptive, - num_neutral_fluids, background_neutral_velocity = background_neutral_velocity(config), background_neutral_density = background_neutral_density(config), Bmax = maximum(cache.B), diff --git a/src/simulation/solution.jl b/src/simulation/solution.jl index 34cc8f0c..af606712 100644 --- a/src/simulation/solution.jl +++ b/src/simulation/solution.jl @@ -35,13 +35,13 @@ end @inline _saved_fields_vector() = ( :μ, :Tev, :ϕ, :∇ϕ, :ne, :pe, :ue, :∇pe, :νan, :νc, :νen, - :νei, :radial_loss_frequency, :νew_momentum, :νiz, :νex, :νe, :Id, :ji, :nn_tot, + :νei, :radial_loss_frequency, :νew_momentum, :νiz, :νex, :νe, :Id, :ji, :nn, :anom_multiplier, :ohmic_heating, :wall_losses, :inelastic_losses, :Vs, :channel_area, :inner_radius, :outer_radius, :dA_dz, :tanδ, :anom_variables, :dt ) -@inline _saved_fields_matrix() = (:ni, :ui, :niui, :nn) +@inline _saved_fields_matrix() = (:ni, :ui, :niui) @inline saved_fields() = (_saved_fields_vector()..., _saved_fields_matrix()...) function solve(U, params, tspan; saveat) diff --git a/src/simulation/sourceterms.jl b/src/simulation/sourceterms.jl index dd3ec0e5..e1373d15 100644 --- a/src/simulation/sourceterms.jl +++ b/src/simulation/sourceterms.jl @@ -33,7 +33,7 @@ function apply_reactions!(dU::AbstractArray{T}, U::AbstractArray{T}, params, i:: if !params.config.LANDMARK # Momentum transfer due to ionization - if reactant_index == index.ρn[1] + if reactant_index == index.ρn reactant_velocity = params.config.neutral_velocity else reactant_velocity = U[reactant_index + 1, i] / U[reactant_index, i] @@ -69,24 +69,17 @@ function apply_ion_acceleration!(dU, U, params, i) end function apply_ion_wall_losses!(dU, U, params, i) - (;index, config, z_cell, L_ch) = params + (;index, config, z_cell, L_ch, cache) = params (;ncharge, propellant, wall_loss_model, thruster) = config geometry = thruster.geometry Δr = geometry.outer_radius - geometry.inner_radius mi = propellant.m - - if wall_loss_model isa WallSheath - α = wall_loss_model.α - elseif wall_loss_model isa NoWallLosses - return - else - α = 1.0 - end + α = wall_loss_coeff(wall_loss_model) @inbounds for Z in 1:ncharge - in_channel = params.config.transition_function(z_cell[i], L_ch, 1.0, 0.0) + in_channel = config.transition_function(z_cell[i], L_ch, 1.0, 0.0) u_bohm = sqrt(Z * e * params.cache.Tev[i] / mi) νiw = α * in_channel * u_bohm / Δr @@ -97,7 +90,7 @@ function apply_ion_wall_losses!(dU, U, params, i) dU[index.ρiui[Z], i] -= momentum_loss # Neutrals gain density due to ion recombination at the walls - dU[index.ρn[1], i] += density_loss + dU[index.ρn, i] += density_loss end end diff --git a/src/simulation/update_electrons.jl b/src/simulation/update_electrons.jl index e10c867a..1fad9ae0 100644 --- a/src/simulation/update_electrons.jl +++ b/src/simulation/update_electrons.jl @@ -3,7 +3,7 @@ function update_electrons!(U, params, t = 0) (;index, control_current, target_current, Kp, Ti, mi) = params (; B, ue, Tev, ∇ϕ, ϕ, pe, ne, nϵ, μ, ∇pe, νan, νc, νen, νei, radial_loss_frequency, - Z_eff, νiz, νex, νe, ji, Id, νew_momentum, κ, ni, ui, Vs, nn, nn_tot, niui, + Z_eff, νiz, νex, νe, ji, Id, νew_momentum, κ, ni, ui, Vs, nn, niui, Id_smoothed, smoothing_time_constant, anom_multiplier, errors, channel_area ) = params.cache @@ -23,11 +23,7 @@ function update_electrons!(U, params, t = 0) # Update plasma quantities @inbounds for i in 1:ncells # Compute number density for each neutral fluid - nn_tot[i] = 0.0 - for j in 1:params.num_neutral_fluids - nn[j, i] = U[index.ρn[j], i] / params.config.propellant.m - nn_tot[i] += nn[j, i] - end + nn[i] = U[index.ρn, i] / params.config.propellant.m # Compute ion derived quantities ne[i] = 0.0 @@ -74,16 +70,16 @@ function update_electrons!(U, params, t = 0) # Update other collisions @inbounds for i in 1:ncells # Compute electron-neutral and electron-ion collision frequencies - νen[i] = freq_electron_neutral(params.electron_neutral_collisions, nn_tot[i], Tev[i]) + νen[i] = freq_electron_neutral(params.electron_neutral_collisions, nn[i], Tev[i]) # Compute total classical collision frequency # If we're not running the LANDMARK benchmark, include momentum transfer due to inelastic collisions νc[i] = νen[i] + νei[i] + !params.config.LANDMARK * (νiz[i] + νex[i]) - # Compute wall collision frequency, with transition function to force no momentum wall collisions in plume + # Compute wall collision frequency, with transition function to force no momentum wall collisions in plume radial_loss_frequency[i] = freq_electron_wall(params.config.wall_loss_model, U, params, i) νew_momentum[i] = radial_loss_frequency[i]* params.config.transition_function(params.z_cell[i], params.L_ch, 1.0, 0.0) - + end # Update anomalous transport diff --git a/src/visualization/recipes.jl b/src/visualization/recipes.jl index 1bcf665b..85f433ef 100644 --- a/src/visualization/recipes.jl +++ b/src/visualization/recipes.jl @@ -36,7 +36,7 @@ # Plot neutral density @series begin - y := sol[:nn_tot][frame] + y := sol[:nn][frame] ylabel := "Density (m⁻³)" subplot := 1 title := "Neutral density" @@ -101,7 +101,7 @@ for rxn in ionization_reactions if rxn.product.Z == Z if rxn.reactant.Z == 0 - reactant_density = sol[:nn_tot][frame] + reactant_density = sol[:nn][frame] else reactant_density = sol[:ni, rxn.reactant.Z][frame] end diff --git a/test/order_verification/ovs_ions.jl b/test/order_verification/ovs_ions.jl index 052d033f..59b30877 100644 --- a/test/order_verification/ovs_ions.jl +++ b/test/order_verification/ovs_ions.jl @@ -136,7 +136,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) z_end = z_cell[end] z_start = z_cell[1] line(v0, v1, z) = v0 + (v1 - v0) * (z - z_start) / (z_end - z_start) - U[index.ρn[1], :] = [line(ρn_func(z_start), ρn_func(z_end), z) for z in z_cell] + U[index.ρn, :] = [line(ρn_func(z_start), ρn_func(z_end), z) for z in z_cell] U[index.ρi[1], :] = [line(ρi_func(z_start), ρi_func(z_end), z) for z in z_cell] U[index.ρiui[1], :] = U[index.ρi[1], :] * ui_func(0.0) @@ -164,7 +164,6 @@ function solve_ions(ncells, scheme; t_end = 1e-4) ionization_reactions, ionization_reactant_indices = [index.ρn], ionization_product_indices = [index.ρi[1]], - num_neutral_fluids = 1, background_neutral_density = 0.0, background_neutral_velocity = 1.0, Δz_cell, Δz_edge, @@ -186,7 +185,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) sol = (;t = [t], u = [U]) - ρn_sim = sol.u[end][index.ρn[1], :] + ρn_sim = sol.u[end][index.ρn, :] ρi_sim = sol.u[end][index.ρi[1], :] ρiui_sim = sol.u[end][index.ρiui[1], :] ui_sim = ρiui_sim ./ ρi_sim diff --git a/test/unit_tests/test_boundary_conditions.jl b/test/unit_tests/test_boundary_conditions.jl index 25d234cd..5c8588e3 100644 --- a/test/unit_tests/test_boundary_conditions.jl +++ b/test/unit_tests/test_boundary_conditions.jl @@ -41,7 +41,6 @@ ϕ = [300.0, 300.0], channel_area = ones(2) * config.thruster.geometry.channel_area ), - num_neutral_fluids = 1, fluids, fluid_ranges, species, @@ -98,13 +97,13 @@ # when ion velocity at left boundary is less than bohm speed, it should be accelerated # to reach bohm speed HallThruster.left_boundary_state!(U_b, U1, params) - @test U_b[index.ρn[1]] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un + @test U_b[index.ρn] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un @test U_b[index.ρiui[1]] / U_b[index.ρi[1]] == -u_bohm_1 @test U_b[index.ρiui[2]] / U_b[index.ρi[2]] == -u_bohm_2 # when ion velocity at left boundary is greater than bohm speed, ions have Neumann BC HallThruster.left_boundary_state!(U_b, U2, params) - @test U_b[index.ρn[1]] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un + @test U_b[index.ρn] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un @test U_b[index.ρiui[1]] == U_2[index.ρiui[1]] @test U_b[index.ρiui[2]] == U_2[index.ρiui[2]] @test U_b[index.ρiui[1]] / U_b[index.ρi[1]] == -2 * u_bohm_1 @@ -112,14 +111,14 @@ # Right boundary condition should be Neumann for all species HallThruster.right_boundary_state!(U_b, U1, params) - @test U_b[index.ρn[1]] ≈ U_1[index.ρn[1]] + @test U_b[index.ρn] ≈ U_1[index.ρn] @test U_b[index.ρi[1]] ≈ U_1[index.ρi[1]] @test U_b[index.ρiui[1]] ≈ U_1[index.ρiui[1]] @test U_b[index.ρi[2]] ≈ U_1[index.ρi[2]] @test U_b[index.ρiui[2]] ≈ U_1[index.ρiui[2]] HallThruster.right_boundary_state!(U_b, U2, params) - @test U_b[index.ρn[1]] ≈ U_2[index.ρn[1]] + @test U_b[index.ρn] ≈ U_2[index.ρn] @test U_b[index.ρi[1]] ≈ U_2[index.ρi[1]] @test U_b[index.ρiui[1]] ≈ U_2[index.ρiui[1]] @test U_b[index.ρi[2]] ≈ U_2[index.ρi[2]] diff --git a/test/unit_tests/test_electrons.jl b/test/unit_tests/test_electrons.jl index ec2522f4..12f89408 100644 --- a/test/unit_tests/test_electrons.jl +++ b/test/unit_tests/test_electrons.jl @@ -22,7 +22,7 @@ mi = HallThruster.Xenon.m index = (ρn = [1], ρi = [2], nϵ = 3) - cache = (;nn_tot = [nn], ne = [ne], B = [B], Tev = [Tev], Z_eff = [1.0], νan = [0.0], κ = [0.0], + cache = (;nn = [nn], ne = [ne], B = [B], Tev = [Tev], Z_eff = [1.0], νan = [0.0], κ = [0.0], μ = μ_e, νc = ν_c, ) c1 = 1/160 diff --git a/test/unit_tests/test_initialization.jl b/test/unit_tests/test_initialization.jl index 2a94a557..fbd50307 100644 --- a/test/unit_tests/test_initialization.jl +++ b/test/unit_tests/test_initialization.jl @@ -39,7 +39,7 @@ HallThruster.initialize!(U, params) - @test abs((U[index.ρn[1], 1] / U[index.ρn[1],end]) - 100) < 1 + @test abs((U[index.ρn, 1] / U[index.ρn,end]) - 100) < 1 ne = [HallThruster.electron_density(U, params, i) for i in 1:ncells+2] ϵ = params.cache.nϵ ./ ne diff --git a/test/unit_tests/test_misc.jl b/test/unit_tests/test_misc.jl index d1598a95..cd98cc6e 100644 --- a/test/unit_tests/test_misc.jl +++ b/test/unit_tests/test_misc.jl @@ -6,7 +6,7 @@ anode_mass_flow_rate = 5u"mg/s", LANDMARK = true, ) - + @test_throws ErrorException HallThruster.run_simulation(Landmark_config; dt=5e-9, duration=4e-9, grid = HallThruster.EvenGrid(2), nsave = 10) config = HallThruster.Config(; @@ -17,7 +17,7 @@ ) @test_logs (:warn, "CFL for adaptive timestepping set higher than stability limit of 0.8. setting CFL to 0.799.") HallThruster.run_simulation(config; dt=5e-9, duration=4e-9, grid = HallThruster.EvenGrid(2), nsave = 10, adaptive = true, CFL = 0.9) - + pressure_config = HallThruster.Config(; thruster = HallThruster.SPT_100, domain = (0.0u"cm", 8.0u"cm"), @@ -25,7 +25,7 @@ anode_mass_flow_rate = 5u"mg/s", background_pressure = 1.0u"Pa", ) - + @test_logs (:warn, "Background neutral pressure set but solve background neutrals not enabled. Did you mean to set solve_background_neutrals to true?") HallThruster.run_simulation(pressure_config; dt=5e-9, duration=4e-9, grid = HallThruster.EvenGrid(2), nsave = 10) end @@ -55,16 +55,16 @@ end end @testset "Stage limiter" begin - index = (ρn = [1], ρi = [2], ρiui = [3], nϵ = 4) + index = (ρn = 1, ρi = [2], ρiui = [3], nϵ = 4) config = (ncharge = 1, min_number_density = 1e6, min_electron_temperature = 1.0, propellant = HallThruster.Xenon) - p = (; config, index, num_neutral_fluids = 1, cache = (;nϵ = [-1.0])) + p = (; config, index, cache = (;nϵ = [-1.0])) U = [-1.0, -1.0, -1.0, -1.0] HallThruster.stage_limiter!(U, p) mi = config.propellant.m - @test U[index.ρn[1]] == config.min_number_density * mi + @test U[index.ρn] == config.min_number_density * mi @test U[index.ρi[1]] == config.min_number_density * mi @test U[index.ρiui[1]] == 1.0 * config.min_number_density * mi @test p.cache.nϵ[1] == config.min_number_density * config.min_electron_temperature @@ -121,13 +121,13 @@ end @test size(U) == (nvars, ncells+2) - (; Aϵ, bϵ, B, νan, νc, μ, ϕ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, F, UL, UR, ni, ui, niui, nn, nn_tot, ji) = cache + (; Aϵ, bϵ, B, νan, νc, μ, ϕ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, F, UL, UR, ni, ui, niui, nn, ji) = cache for arr in (F, UL, UR) @test size(arr) == (nvars, ncells+1) end - for arr in (bϵ, B, νan, νc, μ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, nn_tot, ji) + for arr in (bϵ, B, νan, νc, μ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, ji) @test size(arr) == (ncells+2,) end diff --git a/test/unit_tests/test_walls.jl b/test/unit_tests/test_walls.jl index 7ce95591..a82287ce 100644 --- a/test/unit_tests/test_walls.jl +++ b/test/unit_tests/test_walls.jl @@ -109,7 +109,7 @@ @test HallThruster.freq_electron_wall(landmark_losses, U, params, 3) * params.config.transition_function(params.z_cell[3], params.L_ch, 1.0, 0.0) == 0.0e7 γ = HallThruster.SEE_yield(BN, Tev, γmax) - νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr + νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr νew = νiw / (1 - γ) params.cache.γ_SEE .= γ @@ -259,7 +259,7 @@ end HallThruster.apply_ion_wall_losses!(dU, U, params_constant_sheath, i) # Neutrals should recombine at walls - @test dU[index.ρn[1], i] ≈ -(dU[index.ρi[1], i] + dU[index.ρi[2], i]) + @test dU[index.ρn, i] ≈ -(dU[index.ρi[1], i] + dU[index.ρi[2], i]) # Rate of ion loss is equal to ni νiw Δz = z_edge[2] - z_edge[1] @@ -309,7 +309,7 @@ end HallThruster.apply_ion_wall_losses!(dU, U, params_wall_sheath, i) # Neutrals should recombine at walls - @test dU[index.ρn[1], i] ≈ -(dU[index.ρi[1], i] + dU[index.ρi[2], i]) + @test dU[index.ρn, i] ≈ -(dU[index.ρi[1], i] + dU[index.ρi[2], i]) # Rate of ion loss is equal to Iiw / e / V_cell @test dU[index.ρi[1], i] ≈ -νiw * U[index.ρi[1], i] From 0db013ef33e3d7f27891e9c1fe241cecb16a0809 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 20:53:35 -0400 Subject: [PATCH 02/14] remove dependency on U from almost all electron funcs --- src/simulation/electronenergy.jl | 16 ++--- src/simulation/potential.jl | 11 ++-- src/simulation/restart.jl | 2 +- src/simulation/sourceterms.jl | 26 ++++---- src/simulation/update_electrons.jl | 15 ++--- .../constant_sheath_potential.jl | 9 +-- src/wall_loss_models/no_wall_losses.jl | 7 +-- src/wall_loss_models/wall_losses.jl | 25 ++++---- src/wall_loss_models/wall_sheath.jl | 22 +++---- test/order_verification/ovs_energy.jl | 33 +++++----- test/runtests.jl | 2 +- test/unit_tests/test_initialization.jl | 4 +- test/unit_tests/test_misc.jl | 4 +- test/unit_tests/test_walls.jl | 60 ++++++++----------- 14 files changed, 109 insertions(+), 127 deletions(-) diff --git a/src/simulation/electronenergy.jl b/src/simulation/electronenergy.jl index fc4d7c76..86a4ec34 100644 --- a/src/simulation/electronenergy.jl +++ b/src/simulation/electronenergy.jl @@ -1,11 +1,11 @@ -function update_electron_energy!(U, params, dt) - (;Δz_cell, Δz_edge, index, config, cache, Te_L, Te_R) = params - (;Aϵ, bϵ, nϵ, ue, ne, Tev, channel_area, dA_dz, κ) = cache +function update_electron_energy!(params, dt) + (;Δz_cell, Δz_edge, config, cache, Te_L, Te_R) = params + (;Aϵ, bϵ, nϵ, ue, ne, Tev, channel_area, dA_dz, κ, ni, niui) = cache implicit = params.config.implicit_energy explicit = 1 - implicit - ncells = size(U, 2) + ncells = length(ne) mi = params.config.propellant.m Aϵ.d[1] = 1.0 @@ -25,9 +25,9 @@ function update_electron_energy!(U, params, dt) bϵ[end] = 1.5 * Te_R * ne[end] @inbounds for i in 2:ncells-1 - Q = source_electron_energy(U, params, i) + Q = source_electron_energy(params, i) # User-provided source term - Q += config.source_energy(U, params, i) + Q += config.source_energy(params, i) neL = ne[i-1] ne0 = ne[i] @@ -71,10 +71,10 @@ function update_electron_energy!(U, params, dt) jd = params.cache.Id[] / channel_area[1] # current densities at sheath edge - ji_sheath_edge = e * sum(Z * U[index.ρiui[Z], 1] for Z in 1:params.config.ncharge) / mi + ji_sheath_edge = e * sum(Z * niui[Z, 1] for Z in 1:params.config.ncharge) je_sheath_edge = jd - ji_sheath_edge - ne_sheath_edge = sum(Z * U[index.ρi[Z], 1] for Z in 1:params.config.ncharge) / mi + ne_sheath_edge = sum(Z * ni[Z, 1] for Z in 1:params.config.ncharge) ue_sheath_edge = -je_sheath_edge / ne_sheath_edge / e FL_factor_L = 0.0 diff --git a/src/simulation/potential.jl b/src/simulation/potential.jl index 933197fe..2f5ba683 100644 --- a/src/simulation/potential.jl +++ b/src/simulation/potential.jl @@ -7,17 +7,14 @@ function solve_potential!(ϕ, params) return ϕ end -function anode_sheath_potential(U, params) - (;config, index) = params - (;ne, channel_area) = params.cache - - mi = config.propellant.m +function anode_sheath_potential(params) + (;config) = params + (;ne, niui, channel_area) = params.cache if params.config.LANDMARK return 0.0 end - # Compute anode sheath potential if config.anode_boundary_condition == :sheath ce = sqrt(8 * e * params.cache.Tev[1] / π / me) @@ -27,7 +24,7 @@ function anode_sheath_potential(U, params) jd = params.cache.Id[] / channel_area[1] # current densities at sheath edge - ji_sheath_edge = e * sum(Z * U[index.ρiui[Z], 1] for Z in 1:params.config.ncharge) / mi + ji_sheath_edge = e * sum(Z * niui[Z, 1] for Z in 1:params.config.ncharge) je_sheath_edge = jd - ji_sheath_edge current_ratio = je_sheath_edge / je_sheath diff --git a/src/simulation/restart.jl b/src/simulation/restart.jl index ac57978f..fb9416de 100644 --- a/src/simulation/restart.jl +++ b/src/simulation/restart.jl @@ -93,7 +93,7 @@ function load_restart(grid, fluids, config, sol::Solution) @. @views cache[field][Z, :] = U[2 * Z + 1, :] / config.propellant.m end elseif field == :nn - @. @views cache[field][1, :] = U[1, :] / config.propellant.m + @. @views cache[field] = U[1, :] / config.propellant.m else itp = LinearInterpolation(z_cell, sv[field]) cache[field] .= itp.(grid.cell_centers) diff --git a/src/simulation/sourceterms.jl b/src/simulation/sourceterms.jl index e1373d15..5e212234 100644 --- a/src/simulation/sourceterms.jl +++ b/src/simulation/sourceterms.jl @@ -69,14 +69,20 @@ function apply_ion_acceleration!(dU, U, params, i) end function apply_ion_wall_losses!(dU, U, params, i) - (;index, config, z_cell, L_ch, cache) = params + (;index, config, z_cell, L_ch) = params (;ncharge, propellant, wall_loss_model, thruster) = config geometry = thruster.geometry Δr = geometry.outer_radius - geometry.inner_radius mi = propellant.m - α = wall_loss_coeff(wall_loss_model) + if wall_loss_model isa WallSheath + α = wall_loss_model.α + elseif wall_loss_model isa NoWallLosses + return + else + α = 1.0 + end @inbounds for Z in 1:ncharge in_channel = config.transition_function(z_cell[i], L_ch, 1.0, 0.0) @@ -94,9 +100,9 @@ function apply_ion_wall_losses!(dU, U, params, i) end end -function excitation_losses!(U, params, i) +function excitation_losses!(params, i) (;excitation_reactions, cache, excitation_reactant_indices) = params - (;νex, Tev, inelastic_losses) = cache + (;νex, Tev, inelastic_losses, nn) = cache mi = params.config.propellant.m ne = cache.ne[i] @@ -114,8 +120,8 @@ function excitation_losses!(U, params, i) νex[i] = 0.0 @inbounds for (reactant_inds, rxn) in zip(excitation_reactant_indices, excitation_reactions) rate_coeff = rxn.rate_coeff(ϵ) - for reactant_index in reactant_inds - n_reactant = U[reactant_index, i] / mi + for _ in reactant_inds + n_reactant = nn[i] ndot = reaction_rate(rate_coeff, ne, n_reactant) inelastic_losses[i] += ndot * (rxn.energy - K) νex[i] += ndot / ne @@ -125,7 +131,7 @@ function excitation_losses!(U, params, i) return inelastic_losses[i] end -function electron_kinetic_energy(U, params, i) +function electron_kinetic_energy(params, i) νe = params.cache.νe[i] B = params.cache.B[i] ue = params.cache.ue[i] @@ -133,7 +139,7 @@ function electron_kinetic_energy(U, params, i) return 0.5 * me * (1 + Ωe^2) * ue^2 / e end -function source_electron_energy(U, params, i) +function source_electron_energy(params, i) ne = params.cache.ne[i] ue = params.cache.ue[i] ∇ϕ = params.cache.∇ϕ[i] @@ -154,10 +160,10 @@ function source_electron_energy(U, params, i) end # add excitation losses to total inelastic losses - excitation_losses!(U, params, i) + excitation_losses!(params, i) # compute wall losses - wall_loss = ne * wall_power_loss(params.config.wall_loss_model, U, params, i) + wall_loss = ne * wall_power_loss(params.config.wall_loss_model, params, i) params.cache.wall_losses[i] = wall_loss params.cache.ohmic_heating[i] = ohmic_heating diff --git a/src/simulation/update_electrons.jl b/src/simulation/update_electrons.jl index 1fad9ae0..4141056f 100644 --- a/src/simulation/update_electrons.jl +++ b/src/simulation/update_electrons.jl @@ -77,7 +77,7 @@ function update_electrons!(U, params, t = 0) νc[i] = νen[i] + νei[i] + !params.config.LANDMARK * (νiz[i] + νex[i]) # Compute wall collision frequency, with transition function to force no momentum wall collisions in plume - radial_loss_frequency[i] = freq_electron_wall(params.config.wall_loss_model, U, params, i) + radial_loss_frequency[i] = freq_electron_wall(params.config.wall_loss_model, params, i) νew_momentum[i] = radial_loss_frequency[i]* params.config.transition_function(params.z_cell[i], params.L_ch, 1.0, 0.0) end @@ -100,10 +100,10 @@ function update_electrons!(U, params, t = 0) end # Compute anode sheath potential - Vs[] = anode_sheath_potential(U, params) + Vs[] = anode_sheath_potential(params) # Compute the discharge current by integrating the momentum equation over the whole domain - Id[] = discharge_current(U, params) + Id[] = _discharge_current(params) # Compute the electron velocity and electron kinetic energy @inbounds for i in 1:ncells @@ -111,7 +111,7 @@ function update_electrons!(U, params, t = 0) ue[i] = (ji[i] - Id[] / channel_area[i]) / e / ne[i] # Kinetic energy in both axial and azimuthal directions is accounted for - params.cache.K[i] = electron_kinetic_energy(U, params, i) + params.cache.K[i] = electron_kinetic_energy(params, i) end # Compute potential gradient and pressure gradient @@ -127,7 +127,7 @@ function update_electrons!(U, params, t = 0) params.config.conductivity_model(κ, params) # Update the electron temperature and pressure - update_electron_energy!(U, params, params.dt[]) + update_electron_energy!(params, params.dt[]) # Update the anomalous collision frequency multiplier to match target # discharge current @@ -156,11 +156,12 @@ function update_electrons!(U, params, t = 0) end -function discharge_current(U::Array, params) +# TODO: differentiate this from the postprocessing one +function _discharge_current(params) (;cache, Δz_edge, ϕ_L, ϕ_R) = params (;∇pe, μ, ne, ji, Vs, channel_area) = cache - ncells = size(U, 2) + ncells = length(ne) int1 = 0.0 int2 = 0.0 diff --git a/src/wall_loss_models/constant_sheath_potential.jl b/src/wall_loss_models/constant_sheath_potential.jl index fcfcf634..7f226aa0 100644 --- a/src/wall_loss_models/constant_sheath_potential.jl +++ b/src/wall_loss_models/constant_sheath_potential.jl @@ -4,13 +4,10 @@ Base.@kwdef struct ConstantSheathPotential <: WallLossModel outer_loss_coeff::Float64 end -function freq_electron_wall(::ConstantSheathPotential, U, params, i) - (;z_cell, L_ch) = params - return 1e7 -end +freq_electron_wall(::ConstantSheathPotential, params, i) = 1e7 -function wall_power_loss(model::ConstantSheathPotential, U, params, i) - (;z_cell, index, L_ch, config, cache) = params +function wall_power_loss(model::ConstantSheathPotential, params, i) + (;z_cell, L_ch, config, cache) = params ne = cache.ne[i] nϵ = cache.nϵ[i] diff --git a/src/wall_loss_models/no_wall_losses.jl b/src/wall_loss_models/no_wall_losses.jl index b7bd9805..beeb0a60 100644 --- a/src/wall_loss_models/no_wall_losses.jl +++ b/src/wall_loss_models/no_wall_losses.jl @@ -1,6 +1,5 @@ struct NoWallLosses <: WallLossModel end -freq_electron_wall(model::NoWallLosses, U, params, i) = 0.0 -wall_electron_current(model::NoWallLosses, U, params, i) = 0.0 -wall_power_loss(model::NoWallLosses, U, params, i) = 0.0 - +freq_electron_wall(model::NoWallLosses, params, i) = 0.0 +wall_electron_current(model::NoWallLosses, params, i) = 0.0 +wall_power_loss(model::NoWallLosses, params, i) = 0.0 diff --git a/src/wall_loss_models/wall_losses.jl b/src/wall_loss_models/wall_losses.jl index d94cf329..590390af 100644 --- a/src/wall_loss_models/wall_losses.jl +++ b/src/wall_loss_models/wall_losses.jl @@ -12,38 +12,35 @@ Abstract type for wall loss models in the electron energy equation. Types includ See `WallMaterial` for material options. Users implementing their own `WallLossModel` will need to implement at least three methods - 1) `freq_electron_wall(model, U, params, i)`: Compute the electron-wall momentum transfer collision frequency in cell `i` - 2) `wall_power_loss(model, U, params, i)`: Compute the electron power lost to the walls + 1) `freq_electron_wall(model, params, i)`: Compute the electron-wall momentum transfer collision frequency in cell `i` + 2) `wall_power_loss(model, params, i)`: Compute the electron power lost to the walls -A third method, `wall_electron_current(model, U, params, i)`, will compute the electron current to the walls in cell `i`. If left unimplemented, +A third method, `wall_electron_current(model, params, i)`, will compute the electron current to the walls in cell `i`. If left unimplemented, it defaults to Ie,w = e ne νew V_cell where V_cell is the cell volume. -A fourth method, `wall_ion_current(model, U, params, i, Z)`, for computing the current of ions of charge Z to the walls in cell `i`, may also be implemented. +A fourth method, `wall_ion_current(model, params, i, Z)`, for computing the current of ions of charge Z to the walls in cell `i`, may also be implemented. If left unimplemented, it will default to computing the current assuming Ie,w = Ii,w. """ abstract type WallLossModel end -function freq_electron_wall(model::WallLossModel, U, params, i) +function freq_electron_wall(model::WallLossModel, params, i) error("freq_electron_wall not implemented for wall loss model of type $(typeof(model)). See documentation for WallLossModel for a list of required methods") end -function wall_power_loss(model::WallLossModel, U, params, i) +function wall_power_loss(model::WallLossModel, params, i) error("wall_power_loss not implemented for wall loss model of type $(typeof(model)). See documentation for WallLossModel for a list of required methods") end -function wall_electron_current(::WallLossModel, U, params, i) +function wall_electron_current(::WallLossModel, params, i) (;Δz_cell, cache, A_ch) = params (;ne, νew_momentum) = cache V_cell = A_ch * Δz_cell[i] return e * νew_momentum[i] * V_cell * ne[i] end -function wall_ion_current(model::WallLossModel, U, params, i, Z) - (;index, cache, config) = params - (;ne) = cache - mi = config.propellant.m - ni_Z = U[index.ρi[Z], i] / mi +function wall_ion_current(model::WallLossModel, params, i, Z) + (;ne, ni) = params.cache - Iew = wall_electron_current(model, U, params, i) - return Z * ni_Z / ne[i] * Iew * (1 - params.cache.γ_SEE[i]) + Iew = wall_electron_current(model, params, i) + return Z * ni[Z, i] / ne[i] * Iew * (1 - params.cache.γ_SEE[i]) end diff --git a/src/wall_loss_models/wall_sheath.jl b/src/wall_loss_models/wall_sheath.jl index 23c25012..624c8153 100644 --- a/src/wall_loss_models/wall_sheath.jl +++ b/src/wall_loss_models/wall_sheath.jl @@ -16,7 +16,7 @@ const SiliconDioxide = WallMaterial(name = "Silicon Dioxide", ϵ_star = 18, σ const BNSiO2 = WallMaterial(name = "Boron Nitride Silicon Dioxide", ϵ_star = 40, σ₀ = 0.54) const SiliconCarbide = WallMaterial(name = "Silicon Carbide", ϵ_star = 43, σ₀ = 0.69) -function wall_electron_temperature(U, params, i) +function wall_electron_temperature(params, i) (;cache, config, z_cell) = params shielded = config.thruster.shielded @@ -49,35 +49,35 @@ Base.@kwdef struct WallSheath <: WallLossModel end end -function freq_electron_wall(model::WallSheath, U, params, i) +function freq_electron_wall(model::WallSheath, params, i) (; index, config, cache) = params (; ncharge) = config mi = config.propellant.m - #compute radii difference + #compute radii difference geometry = config.thruster.geometry Δr = geometry.outer_radius - geometry.inner_radius #compute electron wall temperature - Tev = wall_electron_temperature(U, params, i) + Tev = wall_electron_temperature(params, i) #calculate and store SEE coefficient γ = SEE_yield(model.material, Tev, params.γ_SEE_max) cache.γ_SEE[i] = γ - #compute the ion current to the walls + #compute the ion current to the walls j_iw = 0.0 for Z in 1:ncharge - niw = U[index.ρi[Z], i] / mi + niw = cache.ni[Z, i] j_iw += model.α * Z * niw * sqrt(Z * e * Tev / mi) end - #compute electron wall collision frequency + #compute electron wall collision frequency νew = j_iw / (Δr * (1 - γ)) / cache.ne[i] - return νew + return νew end -function wall_power_loss(model::WallSheath, U, params, i) +function wall_power_loss(model::WallSheath, params, i) (;config) = params mi = config.propellant.m - Tev = wall_electron_temperature(U, params, i) + Tev = wall_electron_temperature(params, i) # space charge limited SEE coefficient γ = params.cache.γ_SEE[i] @@ -85,7 +85,7 @@ function wall_power_loss(model::WallSheath, U, params, i) # Space charge-limited sheath potential ϕ_s = sheath_potential(Tev, γ, mi) - # Compute electron wall collision frequency with transition function for energy wall collisions in plume + # Compute electron wall collision frequency with transition function for energy wall collisions in plume νew = params.cache.radial_loss_frequency[i] * params.config.transition_function(params.z_cell[i], params.L_ch, 1.0, params.config.electron_plume_loss_scale) # Compute wall power loss rate diff --git a/test/order_verification/ovs_energy.jl b/test/order_verification/ovs_energy.jl index f821c5e0..42ad5c91 100644 --- a/test/order_verification/ovs_energy.jl +++ b/test/order_verification/ovs_energy.jl @@ -18,8 +18,7 @@ ui = sin_wave(x/L, amplitude = 13000, phase = π/4, nwaves = 0.75, offset = 1000 μ = sin_wave(x/L, amplitude = 1e4, phase = π/2, nwaves = 1.2, offset = 1.1e4) ϵ = sin_wave(x/L, amplitude = 20, phase = 1.3*π/2, nwaves = 1.1, offset = 30) ∇ϕ = Dx(ϕ) -ρiui = ne * ui * HallThruster.Xenon.m -ρn = nn * HallThruster.Xenon.m +niui = ne * ui nϵ = ne * ϵ ue = μ * (∇ϕ - Dx(nϵ)/ne) κ = 10/9 * μ * nϵ @@ -27,12 +26,12 @@ ue = μ * (∇ϕ - Dx(nϵ)/ne) ϕ_func = eval(build_function(ϕ, [x])) ne_func = eval(build_function(ne, [x])) μ_func = eval(build_function(μ, [x])) -ρiui_func = eval(build_function(ρiui, [x])) +niui_func = eval(build_function(niui, [x])) nϵ_func = eval(build_function(nϵ, [x])) κ_func = eval(build_function(κ, [x])) ue_func = eval(build_function(expand_derivatives(ue), [x])) ∇ϕ_func = eval(build_function(expand_derivatives(∇ϕ), [x])) -ρn_func = eval(build_function(ρn, [x])) +nn_func = eval(build_function(nn, [x])) ϵ_func = eval(build_function(ϵ, [x])) k(ϵ) = 8.32 * OVS_rate_coeff_ex(ϵ) @@ -40,14 +39,14 @@ W(ϵ) = 1e7 * ϵ * exp(-20 / ϵ) energy_eq = Dt(nϵ) + Dx(5/3 * nϵ * ue - 10/9 * μ * nϵ * Dx(nϵ/ne)) + ne * (-ue * Dx(ϕ) + nn * k(ϵ) + W(ϵ)) source_energy = eval(build_function(expand_derivatives(energy_eq), [x])) -function solve_energy!(U, params, max_steps, dt, rtol = sqrt(eps(Float64))) +function solve_energy!(params, max_steps, dt, rtol = sqrt(eps(Float64))) t = 0.0 nϵ_old = copy(params.cache.nϵ) residual = Inf iter = 0 res0 = 0.0 while iter < max_steps && abs(residual / res0) > rtol - HallThruster.update_electron_energy!(U, params, dt) + HallThruster.update_electron_energy!(params, dt) params.cache.νiz .= 0.0 params.cache.νex .= 0.0 params.cache.inelastic_losses .= 0.0 @@ -61,12 +60,10 @@ function solve_energy!(U, params, max_steps, dt, rtol = sqrt(eps(Float64))) iter += 1 end - return U, params + return params end function verify_energy(ncells; niters = 20000) - index = (; ρn = 1, nϵ = 2) - grid = HallThruster.generate_grid(HallThruster.SPT_100.geometry, (0.0, 0.05), UnevenGrid(ncells)) z_cell = grid.cell_centers @@ -79,8 +76,8 @@ function verify_energy(ncells; niters = 20000) ϕ = zeros(ncells) ue = ue_func.(z_cell) ∇ϕ = ∇ϕ_func.(z_cell) - ρn = ρn_func.(z_cell) - U = zeros(2, ncells) + nn = nn_func.(z_cell) + niui = niui_func.(z_cell) Tev = ϵ_func.(z_cell) * 2/3 νex = zeros(ncells) νiz = zeros(ncells) @@ -93,7 +90,6 @@ function verify_energy(ncells; niters = 20000) Te_L = nϵ_exact[1] / ne[1] Te_R = nϵ_exact[end] / ne[end] - U[1, :] = ρn nϵ = Te_L * ne # Set initial temp to 3 eV Aϵ = Tridiagonal(ones(ncells-1), ones(ncells), ones(ncells-1)) @@ -101,7 +97,7 @@ function verify_energy(ncells; niters = 20000) min_electron_temperature = 0.1 * min(Te_L, Te_R) - source_func = (U, params, i) -> source_energy(params.z_cell[i]) + source_func = (params, i) -> source_energy(params.z_cell[i]) # Test backward difference implicit solve dt = 1e-6 @@ -142,13 +138,13 @@ function verify_energy(ncells; niters = 20000) cache = (; Aϵ, bϵ, μ, ϕ, ne, ue, ∇ϕ, Tev, pe, νex, νiz, wall_losses, ohmic_heating, inelastic_losses, - channel_area, dA_dz, κ, nϵ + channel_area, dA_dz, κ, nϵ, nn, ni = ne, niui ) Δz_cell, Δz_edge = HallThruster.grid_spacing(grid) params = (; - z_cell, z_edge, index, Te_L = 2/3 * Te_L, Te_R = 2/3 * Te_R, cache, config, + z_cell, z_edge, Te_L = 2/3 * Te_L, Te_R = 2/3 * Te_R, cache, config, dt, L_ch, propellant, ionization_reactions, ionization_reactant_indices, @@ -158,11 +154,10 @@ function verify_energy(ncells; niters = 20000) Δz_cell, Δz_edge ) - solve_energy!(U, params, niters, dt) + solve_energy!(params, niters, dt) results_implicit = (;z = z_cell, exact = nϵ_exact, sim = params.cache.nϵ[:]) # Test crank-nicholson implicit solve - U[1, :] = ρn params.cache.nϵ .= ne * Te_L config = (; @@ -175,7 +170,7 @@ function verify_energy(ncells; niters = 20000) dt = 8 / maximum(abs.(ue)) * (z_cell[2] - z_cell[1]) params = (; - z_cell, z_edge, index, Te_L = 2/3 * Te_L, Te_R = 2/3 * Te_R, cache, config, + z_cell, z_edge, Te_L = 2/3 * Te_L, Te_R = 2/3 * Te_R, cache, config, dt, L_ch, propellant, ionization_reactions, ionization_reactant_indices, @@ -185,7 +180,7 @@ function verify_energy(ncells; niters = 20000) Δz_cell, Δz_edge ) - solve_energy!(U, params, niters, dt) + solve_energy!(params, niters, dt) results_crank_nicholson = (;z = z_cell, exact = nϵ_exact, sim = params.cache.nϵ[:]) return (results_implicit, results_crank_nicholson) diff --git a/test/runtests.jl b/test/runtests.jl index 83b6f906..9e377701 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,7 +27,6 @@ include("unit_tests/test_json.jl") using Symbolics include("order_verification/ovs_funcs.jl") include("order_verification/ovs_energy.jl") -include("order_verification/ovs_ions.jl") @testset "Order verification (electron energy)" begin vfunc = x -> OVS_Energy.verify_energy(x) @@ -49,6 +48,7 @@ include("order_verification/ovs_ions.jl") end end +include("order_verification/ovs_ions.jl") @testset "Order verification (neutrals and ions)" begin refinements = refines(6, 10, 2) diff --git a/test/unit_tests/test_initialization.jl b/test/unit_tests/test_initialization.jl index fbd50307..337341c5 100644 --- a/test/unit_tests/test_initialization.jl +++ b/test/unit_tests/test_initialization.jl @@ -129,7 +129,7 @@ end index = HallThruster.configure_index(fluids, fluid_ranges) @test keys(index) == (:ρn, :ρi, :ρiui) - @test values(index) == ([1], [2, 4, 6], [3, 5, 7]) + @test values(index) == (1, [2, 4, 6], [3, 5, 7]) # load collisions and reactions ionization_reactions = HallThruster._load_reactions(config.ionization_model, unique(species)) @@ -174,7 +174,7 @@ end index = HallThruster.configure_index(fluids, fluid_ranges) @test keys(index) == (:ρn, :ρi, :ρiui) - @test values(index) == ([1], [2, 4, 6], [3, 5, 7]) + @test values(index) == (1, [2, 4, 6], [3, 5, 7]) # load collisions and reactions ionization_reactions = HallThruster._load_reactions(config.ionization_model, unique(species)) diff --git a/test/unit_tests/test_misc.jl b/test/unit_tests/test_misc.jl index cd98cc6e..370adda2 100644 --- a/test/unit_tests/test_misc.jl +++ b/test/unit_tests/test_misc.jl @@ -127,7 +127,7 @@ end @test size(arr) == (nvars, ncells+1) end - for arr in (bϵ, B, νan, νc, μ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, ji) + for arr in (bϵ, B, νan, νc, μ, ∇ϕ, ne, Tev, pe, ue, ∇pe, νen, νei, radial_loss_frequency, νew_momentum, ji, nn) @test size(arr) == (ncells+2,) end @@ -135,8 +135,6 @@ end @test size(arr) == (3, ncells+2) end - @test size(nn) == (2, ncells+2) - @test size(Aϵ) == (ncells+2, ncells+2) @test Aϵ isa SparseArrays.Tridiagonal diff --git a/test/unit_tests/test_walls.jl b/test/unit_tests/test_walls.jl index a82287ce..62fe5688 100644 --- a/test/unit_tests/test_walls.jl +++ b/test/unit_tests/test_walls.jl @@ -17,10 +17,10 @@ nϵ = 3/2 * ne * Tev cache = (; ne = [ne, ne, ne, ne], + ni = [ne ne ne ne], Tev = [Tev, Tev, Tev, Tev], nϵ = [nϵ, nϵ, nϵ, nϵ], Z_eff = [1.0, 1.0, 1.0, 1.0], - ni = [ne ne ne ne], γ_SEE = [0.0, 0.0, 0.0, 0.0], radial_loss_frequency = [0.0, 0.0, 0.0, 0.0], νew_momentum = [0.0, 0.0, 0.0, 0.0], @@ -66,16 +66,9 @@ Δz_edge, Δz_cell, γ_SEE_max = γmax ) - ρi = ne * mi - - - U = [ - ρi ρi ρi ρi - ] - - @test HallThruster.wall_power_loss(no_losses, U, params, 2) == 0.0 - @test HallThruster.wall_power_loss(landmark_losses, U, params, 2) ≈ αin * 6.0 * 1e7 * exp(-20 / 6.0) - @test HallThruster.wall_power_loss(landmark_losses, U, params, 3) ≈ αout * 6.0 * 1e7 * exp(-20 / 6.0) + @test HallThruster.wall_power_loss(no_losses, params, 2) == 0.0 + @test HallThruster.wall_power_loss(landmark_losses, params, 2) ≈ αin * 6.0 * 1e7 * exp(-20 / 6.0) + @test HallThruster.wall_power_loss(landmark_losses, params, 3) ≈ αout * 6.0 * 1e7 * exp(-20 / 6.0) γ1 = 0.5 γ2 = 1.0 @@ -89,7 +82,6 @@ # Check that SEE properly obeys space charge limit at high electron temps @test HallThruster.SEE_yield(BN, 9000.0, γmax) ≈ γmax - @test HallThruster.SEE_yield(BN, 9000.0, γmax_kr) ≈ γmax_kr γ = HallThruster.SEE_yield(BN, Tev, γmax) @@ -102,11 +94,11 @@ Δz = params.z_edge[2] - params.z_edge[1] V_cell = A_ch * Δz - @test HallThruster.freq_electron_wall(no_losses, U, params, 2) == 0.0 - @test HallThruster.freq_electron_wall(no_losses, U, params, 3) == 0.0 + @test HallThruster.freq_electron_wall(no_losses, params, 2) == 0.0 + @test HallThruster.freq_electron_wall(no_losses, params, 3) == 0.0 - @test HallThruster.freq_electron_wall(landmark_losses, U, params, 2) == 1e7 - @test HallThruster.freq_electron_wall(landmark_losses, U, params, 3) * params.config.transition_function(params.z_cell[3], params.L_ch, 1.0, 0.0) == 0.0e7 + @test HallThruster.freq_electron_wall(landmark_losses, params, 2) == 1e7 + @test HallThruster.freq_electron_wall(landmark_losses, params, 3) * params.config.transition_function(params.z_cell[3], params.L_ch, 1.0, 0.0) == 0.0e7 γ = HallThruster.SEE_yield(BN, Tev, γmax) νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr @@ -120,17 +112,17 @@ params.cache.radial_loss_frequency[1:4] .= νew - Iiw = HallThruster.wall_ion_current(sheath_model, U, params, 2, 1) - Iew = HallThruster.wall_electron_current(sheath_model, U, params, 2) + Iiw = HallThruster.wall_ion_current(sheath_model, params, 2, 1) + Iew = HallThruster.wall_electron_current(sheath_model, params, 2) @test Iew ≈ Iiw / (1 - γ) @test Iiw ≈ νiw * HallThruster.e * V_cell * ne @test Iew ≈ νew * HallThruster.e * V_cell * ne - @test HallThruster.freq_electron_wall(sheath_model, U, params, 2) * params.config.transition_function(params.z_cell[2], L_ch, 1.0, 0.0) ≈ νew - @test HallThruster.freq_electron_wall(sheath_model, U, params, 3) * params.config.transition_function(params.z_cell[3], L_ch, 1.0, 0.0) ≈ 0.0 + @test HallThruster.freq_electron_wall(sheath_model, params, 2) * params.config.transition_function(params.z_cell[2], L_ch, 1.0, 0.0) ≈ νew + @test HallThruster.freq_electron_wall(sheath_model, params, 3) * params.config.transition_function(params.z_cell[3], L_ch, 1.0, 0.0) ≈ 0.0 - @test HallThruster.wall_power_loss(sheath_model, U, params, 2) ≈ νew * (2 * (1 - 0.5 * BN.σ₀) * Tev + (1 - γ) * Vs)/γ - @test HallThruster.wall_power_loss(sheath_model, U, params, 4) ≈ 0.0 + @test HallThruster.wall_power_loss(sheath_model, params, 2) ≈ νew * (2 * (1 - 0.5 * BN.σ₀) * Tev + (1 - γ) * Vs)/γ + @test HallThruster.wall_power_loss(sheath_model, params, 4) ≈ 0.0 end @testset "Ion wall losses" begin @@ -247,9 +239,9 @@ end params_constant_sheath = (;base_params..., config = config_constant_sheath) i = 2 - Iiw_1 = HallThruster.wall_ion_current(constant_sheath, U, params_constant_sheath, i, 1) - Iiw_2 = HallThruster.wall_ion_current(constant_sheath, U, params_constant_sheath, i, 2) - Iew = HallThruster.wall_electron_current(constant_sheath, U, params_constant_sheath, i) + Iiw_1 = HallThruster.wall_ion_current(constant_sheath, params_constant_sheath, i, 1) + Iiw_2 = HallThruster.wall_ion_current(constant_sheath, params_constant_sheath, i, 2) + Iew = HallThruster.wall_electron_current(constant_sheath, params_constant_sheath, i) # Ion and electron wall currents are equivalent inside of channel @test Iiw_1 + Iiw_2 == Iew @@ -277,9 +269,9 @@ end # No wall current outside of channel dU .= 0.0 i = 3 - @test HallThruster.wall_ion_current(constant_sheath, U, params_constant_sheath, i, 1) == 0.0 - @test HallThruster.wall_ion_current(constant_sheath, U, params_constant_sheath, i, 2) == 0.0 - @test HallThruster.wall_electron_current(constant_sheath, U, params_constant_sheath, i) == 0.0 + @test HallThruster.wall_ion_current(constant_sheath, params_constant_sheath, i, 1) == 0.0 + @test HallThruster.wall_ion_current(constant_sheath, params_constant_sheath, i, 2) == 0.0 + @test HallThruster.wall_electron_current(constant_sheath, params_constant_sheath, i) == 0.0 HallThruster.apply_ion_wall_losses!(dU, U, params_constant_sheath, 3) @test all(dU[:, 3] .≈ 0.0) @@ -298,9 +290,9 @@ end params_wall_sheath.cache.νew_momentum[1:2] .= νew i = 2 params_wall_sheath.cache.γ_SEE[i] = γ - Iiw_1 = HallThruster.wall_ion_current(wall_sheath, U, params_wall_sheath, i, 1) - Iiw_2 = HallThruster.wall_ion_current(wall_sheath, U, params_wall_sheath, i, 2) - Iew = HallThruster.wall_electron_current(wall_sheath, U, params_wall_sheath, i) + Iiw_1 = HallThruster.wall_ion_current(wall_sheath, params_wall_sheath, i, 1) + Iiw_2 = HallThruster.wall_ion_current(wall_sheath, params_wall_sheath, i, 2) + Iew = HallThruster.wall_electron_current(wall_sheath, params_wall_sheath, i) @test Iiw_1 ≈ Iew * ni_1 / ne * (1 - γ) @test Iiw_2 ≈ 2 * Iew * ni_2 / ne * (1 - γ) @@ -321,9 +313,9 @@ end # No wall losses in plume i = 3 - @test HallThruster.wall_ion_current(wall_sheath, U, params_wall_sheath, i, 1) == 0.0 - @test HallThruster.wall_ion_current(wall_sheath, U, params_wall_sheath, i, 2) == 0.0 - @test HallThruster.wall_electron_current(wall_sheath, U, params_wall_sheath, i) == 0.0 + @test HallThruster.wall_ion_current(wall_sheath, params_wall_sheath, i, 1) == 0.0 + @test HallThruster.wall_ion_current(wall_sheath, params_wall_sheath, i, 2) == 0.0 + @test HallThruster.wall_electron_current(wall_sheath, params_wall_sheath, i) == 0.0 HallThruster.apply_ion_wall_losses!(dU, U, params_wall_sheath, 3) @test all(dU[:, 3] .≈ 0.0) From 5e9d1d3faf43e5f9e8da181997686a3c70dfb08f Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 21:01:58 -0400 Subject: [PATCH 03/14] add ncells to params --- src/numerics/flux.jl | 6 ++---- src/numerics/limiters.jl | 22 +++++++++++----------- src/simulation/electronenergy.jl | 4 +--- src/simulation/heavy_species.jl | 4 +--- src/simulation/plume.jl | 7 ++----- src/simulation/simulation.jl | 1 + src/simulation/update_electrons.jl | 12 +++--------- 7 files changed, 21 insertions(+), 35 deletions(-) diff --git a/src/numerics/flux.jl b/src/numerics/flux.jl index 860ba879..e927bd03 100644 --- a/src/numerics/flux.jl +++ b/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, 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 diff --git a/src/numerics/limiters.jl b/src/numerics/limiters.jl index 2e935b05..364fb5cf 100644 --- a/src/numerics/limiters.jl +++ b/src/numerics/limiters.jl @@ -17,19 +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 - U[p.index.ρn, i] = max(U[p.index.ρn, i], min_density) + 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 diff --git a/src/simulation/electronenergy.jl b/src/simulation/electronenergy.jl index 86a4ec34..c5debf08 100644 --- a/src/simulation/electronenergy.jl +++ b/src/simulation/electronenergy.jl @@ -1,12 +1,10 @@ function update_electron_energy!(params, dt) - (;Δz_cell, Δz_edge, config, cache, Te_L, Te_R) = params + (;Δz_cell, Δz_edge, config, cache, Te_L, Te_R, ncells) = params (;Aϵ, bϵ, nϵ, ue, ne, Tev, channel_area, dA_dz, κ, ni, niui) = cache implicit = params.config.implicit_energy explicit = 1 - implicit - ncells = length(ne) - mi = params.config.propellant.m Aϵ.d[1] = 1.0 Aϵ.du[1] = 0.0 diff --git a/src/simulation/heavy_species.jl b/src/simulation/heavy_species.jl index 5725db23..b8540361 100644 --- a/src/simulation/heavy_species.jl +++ b/src/simulation/heavy_species.jl @@ -1,7 +1,7 @@ function update_heavy_species!(dU, U, params, t; apply_boundary_conditions = true) - (;index, Δz_cell, config, cache) = params + (;index, Δz_cell, config, cache, ncells) = params (; source_neutrals, source_ion_continuity, source_ion_momentum, ncharge, ion_wall_losses @@ -9,8 +9,6 @@ function update_heavy_species!(dU, U, params, t; apply_boundary_conditions = tru (;F, UL, UR, channel_area, dA_dz) = cache - ncells = size(U, 2) - compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions) params.cache.dt[] = Inf diff --git a/src/simulation/plume.jl b/src/simulation/plume.jl index e50fde20..ce4b373f 100644 --- a/src/simulation/plume.jl +++ b/src/simulation/plume.jl @@ -1,14 +1,11 @@ function update_plume_geometry!(U, params; initialize = false) - - (; mi, z_cell, L_ch, exit_plane_index, config, A_ch, cache, index) = params + (; mi, z_cell, L_ch, exit_plane_index, config, A_ch, cache, index, ncells) = params (; channel_area, inner_radius, outer_radius, channel_height, dA_dz, tanδ, Tev) = cache - if !initialize && !params.config.solve_plume return end - ncells = length(z_cell) geometry = config.thruster.geometry r_in = geometry.inner_radius r_out = geometry.outer_radius @@ -39,4 +36,4 @@ function update_plume_geometry!(U, params; initialize = false) end return nothing -end \ No newline at end of file +end diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 39685a1f..f7c994e9 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -205,6 +205,7 @@ function run_simulation( # Simulation parameters params = (; + ncells = ncells+2, ncharge = config.ncharge, mi, config = config, diff --git a/src/simulation/update_electrons.jl b/src/simulation/update_electrons.jl index 4141056f..41d8be16 100644 --- a/src/simulation/update_electrons.jl +++ b/src/simulation/update_electrons.jl @@ -1,6 +1,6 @@ # update useful quantities relevant for potential, electron energy and fluid solve function update_electrons!(U, params, t = 0) - (;index, control_current, target_current, Kp, Ti, mi) = params + (;index, control_current, target_current, Kp, Ti, mi, ncells) = params (; B, ue, Tev, ∇ϕ, ϕ, pe, ne, nϵ, μ, ∇pe, νan, νc, νen, νei, radial_loss_frequency, Z_eff, νiz, νex, νe, ji, Id, νew_momentum, κ, ni, ui, Vs, nn, niui, @@ -18,8 +18,6 @@ function update_electrons!(U, params, t = 0) @views left_boundary_state!(U[:, 1], U, params) @views right_boundary_state!(U[:, end], U, params) - ncells = size(U, 2) - # Update plasma quantities @inbounds for i in 1:ncells # Compute number density for each neutral fluid @@ -158,11 +156,9 @@ end # TODO: differentiate this from the postprocessing one function _discharge_current(params) - (;cache, Δz_edge, ϕ_L, ϕ_R) = params + (;cache, Δz_edge, ϕ_L, ϕ_R, ncells) = params (;∇pe, μ, ne, ji, Vs, channel_area) = cache - ncells = length(ne) - int1 = 0.0 int2 = 0.0 @@ -200,9 +196,7 @@ end function compute_pressure_gradient!(∇pe, params) (; pe) = params.cache - (;z_cell) = params - - ncells = length(z_cell) + (;z_cell, ncells) = params # Pressure gradient (forward) ∇pe[1] = forward_difference(pe[1], pe[2], pe[3], z_cell[1], z_cell[2], z_cell[3]) From 4dd28b08d75539566124bc913bfe3580eb11c065 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 21:25:21 -0400 Subject: [PATCH 04/14] update_electrons no longer relies on U, add update_ions --- src/simulation/heavy_species.jl | 2 +- src/simulation/simulation.jl | 7 ++- src/simulation/solution.jl | 88 +++++++++++++++++++++++------- src/simulation/update_electrons.jl | 53 ++---------------- 4 files changed, 80 insertions(+), 70 deletions(-) diff --git a/src/simulation/heavy_species.jl b/src/simulation/heavy_species.jl index b8540361..b620b3e7 100644 --- a/src/simulation/heavy_species.jl +++ b/src/simulation/heavy_species.jl @@ -1,5 +1,5 @@ -function update_heavy_species!(dU, U, params, t; apply_boundary_conditions = true) +function update_heavy_species!(dU, U, params; apply_boundary_conditions = true) (;index, Δz_cell, config, cache, ncells) = params (; diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index f7c994e9..722f3b9d 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -235,7 +235,9 @@ function run_simulation( Δz_cell, Δz_edge, control_current, target_current, Kp, Ti, Td, exit_plane_index = findfirst(>=(config.thruster.geometry.channel_length), z_cell) - 1, - dtmin, dtmax + dtmin, dtmax, + # landmark benchmark uses pe = 3/2 ne Te, otherwise use pe = ne Te + pe_factor = config.LANDMARK ? 3/2 : 1.0 ) # Compute maximum allowed iterations @@ -250,7 +252,8 @@ function run_simulation( end # make values in params available for first timestep - update_electrons!(U, params) + update_ions!(U, params) + update_electrons!(params) # Set up sol_info = @timed solve(U, params, tspan; saveat) diff --git a/src/simulation/solution.jl b/src/simulation/solution.jl index af606712..2132e098 100644 --- a/src/simulation/solution.jl +++ b/src/simulation/solution.jl @@ -1,20 +1,3 @@ -# Perform one step of the Strong-stability-preserving RK22 algorithm -function ssprk22_step!(U, f, params, t, dt) - (;k, u1) = params.cache - - # First step of SSPRK22 - f(k, U, params, t) - @. u1 = U + dt * k - stage_limiter!(u1, params) - - # Second step of SSPRK22 - f(k, u1, params, t+dt) - @. U = (U + u1 + dt * k) / 2 - stage_limiter!(U, params) - - return nothing -end - struct Solution{T, U, P, S} t::T u::U @@ -33,6 +16,65 @@ function Base.show(io::IO, mime::MIME"text/plain", sol::Solution) print(io, "End time: $(sol.t[end]) seconds") end +# Perform one step of the Strong-stability-preserving RK22 algorithm +function ssprk22_step!(U, f, params, dt) + (;k, u1) = params.cache + + # First step of SSPRK22 + f(k, U, params) + @. u1 = U + dt * k + stage_limiter!(u1, params) + + # Second step of SSPRK22 + f(k, u1, params) + @. U = (U + u1 + dt * k) / 2 + stage_limiter!(U, params) + + return nothing +end + + +function update_ions!(U, params) + (;index, ncells, cache) = params + (;nn, ne, ni, ui, niui, Z_eff, ji) = cache + mi = params.config.propellant.m + + # Update heavy species + ssprk22_step!(U, update_heavy_species!, params, params.dt[]) + + # Apply fluid boundary conditions + @views left_boundary_state!(U[:, 1], U, params) + @views right_boundary_state!(U[:, end], U, params) + + # Update plasma quantities + @inbounds for i in 1:ncells + # Compute number density for each neutral fluid + nn[i] = U[index.ρn, i] / params.config.propellant.m + + # Compute ion derived quantities + ne[i] = 0.0 + Z_eff[i] = 0.0 + ji[i] = 0.0 + @inbounds for Z in 1:params.config.ncharge + _ni = U[index.ρi[Z], i] / mi + _niui = U[index.ρiui[Z], i] / mi + ni[Z, i] = _ni + niui[Z, i] = _niui + ui[Z, i] = _niui / _ni + ne[i] += Z * _ni + Z_eff[i] += _ni + ji[i] += Z * e * _niui + end + + # Compute electron number density, making sure it is above floor + ne[i] = max(params.config.min_number_density, ne[i]) + + # Effective ion charge state (density-weighted average charge state) + Z_eff[i] = max(1.0, ne[i] / Z_eff[i]) + end +end + + @inline _saved_fields_vector() = ( :μ, :Tev, :ϕ, :∇ϕ, :ne, :pe, :ue, :∇pe, :νan, :νc, :νen, :νei, :radial_loss_frequency, :νew_momentum, :νiz, :νex, :νe, :Id, :ji, :nn, @@ -67,15 +109,21 @@ function solve(U, params, tspan; saveat) t += params.dt[] - # Update heavy species - ssprk22_step!(U, update_heavy_species!, params, t, params.dt[]) + # update heavy species quantities + update_ions!(U, params) # Update electron quantities - update_electrons!(U, params, t) + update_electrons!(params, t) # Update plume geometry update_plume_geometry!(U, params) + # Allow for system interrupts + yield() + + # Update the current iteration + params.iteration[1] += 1 + # Check for NaNs and terminate if necessary nandetected = false infdetected = false diff --git a/src/simulation/update_electrons.jl b/src/simulation/update_electrons.jl index 41d8be16..ffa0a244 100644 --- a/src/simulation/update_electrons.jl +++ b/src/simulation/update_electrons.jl @@ -1,61 +1,20 @@ + # update useful quantities relevant for potential, electron energy and fluid solve -function update_electrons!(U, params, t = 0) - (;index, control_current, target_current, Kp, Ti, mi, ncells) = params +function update_electrons!(params, t = 0) + (;control_current, target_current, Kp, Ti, pe_factor, ncells) = params (; B, ue, Tev, ∇ϕ, ϕ, pe, ne, nϵ, μ, ∇pe, νan, νc, νen, νei, radial_loss_frequency, - Z_eff, νiz, νex, νe, ji, Id, νew_momentum, κ, ni, ui, Vs, nn, niui, + Z_eff, νiz, νex, νe, ji, Id, νew_momentum, κ, Vs, nn, Id_smoothed, smoothing_time_constant, anom_multiplier, errors, channel_area ) = params.cache - # Allow for system interrupts - yield() - - # Update the current iteration - params.iteration[1] += 1 - - # Apply fluid boundary conditions - @views left_boundary_state!(U[:, 1], U, params) - @views right_boundary_state!(U[:, end], U, params) - - # Update plasma quantities + # Update plasma quantities based on new density @inbounds for i in 1:ncells - # Compute number density for each neutral fluid - nn[i] = U[index.ρn, i] / params.config.propellant.m - - # Compute ion derived quantities - ne[i] = 0.0 - Z_eff[i] = 0.0 - ji[i] = 0.0 - @inbounds for Z in 1:params.config.ncharge - _ni = U[index.ρi[Z], i] / mi - _niui = U[index.ρiui[Z], i] / mi - ni[Z, i] = _ni - niui[Z, i] = _niui - ui[Z, i] = _niui / _ni - ne[i] += Z * _ni - Z_eff[i] += _ni - ji[i] += Z * e * _niui - end - - # Compute electron number density, making sure it is above floor - ne[i] = max(params.config.min_number_density, ne[i]) - - # Effective ion charge state (density-weighted average charge state) - Z_eff[i] = max(1.0, ne[i] / Z_eff[i]) - # Compute new electron temperature Tev[i] = 2/3 * max(params.config.min_electron_temperature, nϵ[i]/ne[i]) - # Compute electron pressure - pe[i] = if params.config.LANDMARK - # The LANDMARK benchmark uses nϵ instead of pe in the potential solver, but we use pe, so - # we need to define pe = 3/2 ne Tev - 3/2 * ne[i] * Tev[i] - else - # Otherwise, just use typical ideal gas law. - ne[i] * Tev[i] - end + pe[i] = pe_factor * ne[i] * Tev[i] end # Update electron-ion collisions From ae1617cc303c9da5c0a5f7f8ca4ae240d7da8a3e Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 21:35:28 -0400 Subject: [PATCH 05/14] rearrange ion stuff --- src/HallThruster.jl | 4 +- src/simulation/simulation.jl | 2 +- src/simulation/solution.jl | 62 +------------------ ...avy_species.jl => update_heavy_species.jl} | 58 ++++++++++++++++- test/order_verification/ovs_ions.jl | 5 +- 5 files changed, 61 insertions(+), 70 deletions(-) rename src/simulation/{heavy_species.jl => update_heavy_species.jl} (54%) diff --git a/src/HallThruster.jl b/src/HallThruster.jl index 0db7f66a..268962c6 100644 --- a/src/HallThruster.jl +++ b/src/HallThruster.jl @@ -55,12 +55,12 @@ 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/restart.jl") diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 722f3b9d..0a11c160 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -252,7 +252,7 @@ function run_simulation( end # make values in params available for first timestep - update_ions!(U, params) + update_heavy_species!(U, params) update_electrons!(params) # Set up diff --git a/src/simulation/solution.jl b/src/simulation/solution.jl index 2132e098..5c4ce125 100644 --- a/src/simulation/solution.jl +++ b/src/simulation/solution.jl @@ -16,65 +16,6 @@ function Base.show(io::IO, mime::MIME"text/plain", sol::Solution) print(io, "End time: $(sol.t[end]) seconds") end -# Perform one step of the Strong-stability-preserving RK22 algorithm -function ssprk22_step!(U, f, params, dt) - (;k, u1) = params.cache - - # First step of SSPRK22 - f(k, U, params) - @. u1 = U + dt * k - stage_limiter!(u1, params) - - # Second step of SSPRK22 - f(k, u1, params) - @. U = (U + u1 + dt * k) / 2 - stage_limiter!(U, params) - - return nothing -end - - -function update_ions!(U, params) - (;index, ncells, cache) = params - (;nn, ne, ni, ui, niui, Z_eff, ji) = cache - mi = params.config.propellant.m - - # Update heavy species - ssprk22_step!(U, update_heavy_species!, params, params.dt[]) - - # Apply fluid boundary conditions - @views left_boundary_state!(U[:, 1], U, params) - @views right_boundary_state!(U[:, end], U, params) - - # Update plasma quantities - @inbounds for i in 1:ncells - # Compute number density for each neutral fluid - nn[i] = U[index.ρn, i] / params.config.propellant.m - - # Compute ion derived quantities - ne[i] = 0.0 - Z_eff[i] = 0.0 - ji[i] = 0.0 - @inbounds for Z in 1:params.config.ncharge - _ni = U[index.ρi[Z], i] / mi - _niui = U[index.ρiui[Z], i] / mi - ni[Z, i] = _ni - niui[Z, i] = _niui - ui[Z, i] = _niui / _ni - ne[i] += Z * _ni - Z_eff[i] += _ni - ji[i] += Z * e * _niui - end - - # Compute electron number density, making sure it is above floor - ne[i] = max(params.config.min_number_density, ne[i]) - - # Effective ion charge state (density-weighted average charge state) - Z_eff[i] = max(1.0, ne[i] / Z_eff[i]) - end -end - - @inline _saved_fields_vector() = ( :μ, :Tev, :ϕ, :∇ϕ, :ne, :pe, :ue, :∇pe, :νan, :νc, :νen, :νei, :radial_loss_frequency, :νew_momentum, :νiz, :νex, :νe, :Id, :ji, :nn, @@ -110,7 +51,8 @@ function solve(U, params, tspan; saveat) t += params.dt[] # update heavy species quantities - update_ions!(U, params) + integrate_heavy_species!(U, params, params.dt[]) + update_heavy_species!(U, params) # Update electron quantities update_electrons!(params, t) diff --git a/src/simulation/heavy_species.jl b/src/simulation/update_heavy_species.jl similarity index 54% rename from src/simulation/heavy_species.jl rename to src/simulation/update_heavy_species.jl index b620b3e7..05b94c31 100644 --- a/src/simulation/heavy_species.jl +++ b/src/simulation/update_heavy_species.jl @@ -1,6 +1,4 @@ - -function update_heavy_species!(dU, U, params; apply_boundary_conditions = true) - +function iterate_heavy_species!(dU, U, params; apply_boundary_conditions = true) (;index, Δz_cell, config, cache, ncells) = params (; source_neutrals, source_ion_continuity, source_ion_momentum, @@ -66,3 +64,57 @@ function update_heavy_species!(dU, U, params; apply_boundary_conditions = true) return nothing end + +# Perform one step of the Strong-stability-preserving RK22 algorithm to the ion fluid +function integrate_heavy_species!(U, params, dt) + (;k, u1) = params.cache + + # First step of SSPRK22 + iterate_heavy_species!(k, U, params) + @. u1 = U + dt * k + stage_limiter!(u1, params) + + # Second step of SSPRK22 + iterate_heavy_species!(k, u1, params) + @. U = (U + u1 + dt * k) / 2 + stage_limiter!(U, params) + + return nothing +end + +function update_heavy_species!(U, params) + (;index, ncells, cache) = params + (;nn, ne, ni, ui, niui, Z_eff, ji) = cache + mi = params.config.propellant.m + + # Apply fluid boundary conditions + @views left_boundary_state!(U[:, 1], U, params) + @views right_boundary_state!(U[:, end], U, params) + + # Update plasma quantities + @inbounds for i in 1:ncells + # Compute number density for each neutral fluid + nn[i] = U[index.ρn, i] / params.config.propellant.m + + # Compute ion derived quantities + ne[i] = 0.0 + Z_eff[i] = 0.0 + ji[i] = 0.0 + @inbounds for Z in 1:params.config.ncharge + _ni = U[index.ρi[Z], i] / mi + _niui = U[index.ρiui[Z], i] / mi + ni[Z, i] = _ni + niui[Z, i] = _niui + ui[Z, i] = _niui / _ni + ne[i] += Z * _ni + Z_eff[i] += _ni + ji[i] += Z * e * _niui + end + + # Compute electron number density, making sure it is above floor + ne[i] = max(params.config.min_number_density, ne[i]) + + # Effective ion charge state (density-weighted average charge state) + Z_eff[i] = max(1.0, ne[i] / Z_eff[i]) + end +end diff --git a/test/order_verification/ovs_ions.jl b/test/order_verification/ovs_ions.jl index 59b30877..b008962e 100644 --- a/test/order_verification/ovs_ions.jl +++ b/test/order_verification/ovs_ions.jl @@ -176,10 +176,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) t = 0.0 while t < t_end @views U[:, end] = U[:, end-1] - HallThruster.ssprk22_step!( - U, (dU, U, params, t) -> HallThruster.update_heavy_species!(dU, U, params, t; apply_boundary_conditions = false), - params, t, dt[] - ) + HallThruster.integrate_heavy_species!(U, params, dt[]); t += dt[] end From 9033a1b66400bbea2e11ce9df2ac1f8ba917fb59 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 21:51:42 -0400 Subject: [PATCH 06/14] work on making tests pass --- reactions/elastic_Kr.dat | 1 + reactions/elastic_Xe.dat | 1 + src/collisions/collision_frequencies.jl | 2 +- src/simulation/simulation.jl | 6 ++--- src/simulation/solution.jl | 32 ++++++++++++------------- test/unit_tests/test_electrons.jl | 6 ++--- test/unit_tests/test_geometry.jl | 15 ++++++------ test/unit_tests/test_misc.jl | 4 ++-- test/unit_tests/test_restarts.jl | 11 ++++----- 9 files changed, 38 insertions(+), 40 deletions(-) diff --git a/reactions/elastic_Kr.dat b/reactions/elastic_Kr.dat index 6c6f6c3c..849bcbfd 100644 --- a/reactions/elastic_Kr.dat +++ b/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 diff --git a/reactions/elastic_Xe.dat b/reactions/elastic_Xe.dat index 8dfb7849..10693cd6 100644 --- a/reactions/elastic_Xe.dat +++ b/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 diff --git a/src/collisions/collision_frequencies.jl b/src/collisions/collision_frequencies.jl index 73c084b5..1d56fa11 100644 --- a/src/collisions/collision_frequencies.jl +++ b/src/collisions/collision_frequencies.jl @@ -10,7 +10,7 @@ Effective frequency of electron scattering caused by collisions with neutrals return νen end -function freq_electron_neutral(U::AbstractArray, params::NamedTuple, i::Int) +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) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 0a11c160..96589c3b 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -128,12 +128,12 @@ function run_simulation( ) #check that Landmark uses the correct thermal conductivity - if config.LANDMARK & (config.conductivity_model != LANDMARK_conductivity()) + if config.LANDMARK && !(config.conductivity_model isa LANDMARK_conductivity) error("LANDMARK configuration needs to use the LANDMARK thermal conductivity model.") end #check that user is aware of neutral backflow flag - if (config.background_pressure > 0.0) & (config.solve_background_neutrals == false) + if (config.background_pressure > 0.0) && (config.solve_background_neutrals == false) @warn("Background neutral pressure set but solve background neutrals not enabled. Did you mean to set solve_background_neutrals to true?") end @@ -194,7 +194,7 @@ function run_simulation( # this limit is mainly due to empirical testing, but there # may be an analytical reason the ionization timestep cannot use a CFL >= 0.8 if CFL >= 0.8 - @warn("CFL for adaptive timestepping set higher than stability limit of 0.8. setting CFL to 0.799.") + @warn("CFL for adaptive timestepping set higher than stability limit of 0.8. Setting CFL to 0.799.") CFL = 0.799 end end diff --git a/src/simulation/solution.jl b/src/simulation/solution.jl index 5c4ce125..46997605 100644 --- a/src/simulation/solution.jl +++ b/src/simulation/solution.jl @@ -54,18 +54,6 @@ function solve(U, params, tspan; saveat) integrate_heavy_species!(U, params, params.dt[]) update_heavy_species!(U, params) - # Update electron quantities - update_electrons!(params, t) - - # Update plume geometry - update_plume_geometry!(U, params) - - # Allow for system interrupts - yield() - - # Update the current iteration - params.iteration[1] += 1 - # Check for NaNs and terminate if necessary nandetected = false infdetected = false @@ -84,6 +72,22 @@ function solve(U, params, tspan; saveat) end end + if nandetected || infdetected + break + end + + # Update electron quantities + update_electrons!(params, t) + + # Update plume geometry + update_plume_geometry!(U, params) + + # Allow for system interrupts + yield() + + # Update the current iteration + params.iteration[1] += 1 + # Save values at designated intervals # TODO interpolate these to be exact and make a bit more elegant if t > saveat[save_ind] @@ -111,10 +115,6 @@ function solve(U, params, tspan; saveat) save_ind += 1 end - - if nandetected || infdetected - break - end end ind = min(save_ind, length(savevals)+1)-1 diff --git a/test/unit_tests/test_electrons.jl b/test/unit_tests/test_electrons.jl index 12f89408..d57dda4c 100644 --- a/test/unit_tests/test_electrons.jl +++ b/test/unit_tests/test_electrons.jl @@ -51,9 +51,9 @@ params_none = (;cache, index, config = config_none, z_cell = [0.03], L_ch = thruster.geometry.channel_length, electron_neutral_collisions = en_none) U = [mi * nn; mi * ne; ne * 3/2 * Tev ;;] - @test HallThruster.freq_electron_neutral(U, params_landmark, 1) == 2.5e-13 * nn - @test HallThruster.freq_electron_neutral(U, params_gk, 1) == HallThruster.σ_en(Tev) * nn * sqrt(8 * e * Tev / π / me) - @test HallThruster.freq_electron_neutral(U, params_none, 1) == 0.0 + @test HallThruster.freq_electron_neutral(params_landmark, 1) == 2.5e-13 * nn + @test HallThruster.freq_electron_neutral(params_gk, 1) == HallThruster.σ_en(Tev) * nn * sqrt(8 * e * Tev / π / me) + @test HallThruster.freq_electron_neutral(params_none, 1) == 0.0 Z = 1 @test HallThruster.freq_electron_ion(ne, Tev, Z) == 2.9e-12 * Z^2 * ne * HallThruster.coulomb_logarithm(ne, Tev, Z) / Tev^1.5 diff --git a/test/unit_tests/test_geometry.jl b/test/unit_tests/test_geometry.jl index c0006230..e891ce7f 100644 --- a/test/unit_tests/test_geometry.jl +++ b/test/unit_tests/test_geometry.jl @@ -25,15 +25,14 @@ @test HallThruster.channel_perimeter(r1, r0) ≈ 2π * (r1 + r0) @test HallThruster.channel_perimeter(HallThruster.SPT_100) ≈ 2π * (r1 + r0) - - U = zeros(6, 4) + U = zeros(6, 4) U[1, :] .= 1e18 U[2, :] .= 1e17 U[3, :] .= 1e16 U[4, :] = 1e18 .* [-100, 100, 5000, 20000] U[5, :] = 1e17 .* [-100, 100, 2 * 5000, 2 * 20000] U[6, :] = 1e16 .* [-100, 100, 3 * 5000, 3 *20000] - + inner_radius = zeros(4) outer_radius = zeros(4) @@ -46,14 +45,14 @@ tan_δ = 0.5 * max(0.0, min(π/4, cs / ui)) A = π * ((r1 + tan_δ * 0.025)^2 - (max(0.0,r0 - tan_δ * 0.025))^2) - cache = (; channel_area = zeros(4), inner_radius = zeros(4), - outer_radius = zeros(4), channel_height = zeros(4), + cache = (; channel_area = zeros(4), inner_radius = zeros(4), + outer_radius = zeros(4), channel_height = zeros(4), dA_dz = zeros(4), tanδ = zeros(4), Tev = [30]) config = (; thruster = SPT_100, ncharge = 3, solve_plume = true) index = (; ρi = [1,2,3], ρiui = [4,5,6]) params = (; mi = HallThruster.Xenon.m, z_cell = [0.0, 0.01, 0.025, 0.05], - L_ch = 0.025, exit_plane_index = 1, config = config, A_ch = A_ch, - cache = cache, index = index ) + L_ch = 0.025, exit_plane_index = 1, config = config, A_ch = A_ch, + cache = cache, index = index, ncells = size(U, 2)) HallThruster.update_plume_geometry!(U, params) @@ -61,4 +60,4 @@ @test params.cache.dA_dz[2] ≈ 0.0 @test params.cache.dA_dz[3] ≈ 0.0 @test params.cache.dA_dz[4] ≈ (A - A_ch) / 0.025 -end \ No newline at end of file +end diff --git a/test/unit_tests/test_misc.jl b/test/unit_tests/test_misc.jl index 370adda2..979a58c8 100644 --- a/test/unit_tests/test_misc.jl +++ b/test/unit_tests/test_misc.jl @@ -16,7 +16,7 @@ anode_mass_flow_rate = 5u"mg/s", ) - @test_logs (:warn, "CFL for adaptive timestepping set higher than stability limit of 0.8. setting CFL to 0.799.") HallThruster.run_simulation(config; dt=5e-9, duration=4e-9, grid = HallThruster.EvenGrid(2), nsave = 10, adaptive = true, CFL = 0.9) + @test_logs (:warn, "CFL for adaptive timestepping set higher than stability limit of 0.8. Setting CFL to 0.799.") HallThruster.run_simulation(config; dt=5e-9, duration=0e-9, grid = HallThruster.EvenGrid(2), nsave = 10, adaptive = true, CFL = 0.9) pressure_config = HallThruster.Config(; thruster = HallThruster.SPT_100, @@ -58,7 +58,7 @@ end index = (ρn = 1, ρi = [2], ρiui = [3], nϵ = 4) config = (ncharge = 1, min_number_density = 1e6, min_electron_temperature = 1.0, propellant = HallThruster.Xenon) - p = (; config, index, cache = (;nϵ = [-1.0])) + p = (; config, index, cache = (;nϵ = [-1.0]), ncells = 1) U = [-1.0, -1.0, -1.0, -1.0] HallThruster.stage_limiter!(U, p) diff --git a/test/unit_tests/test_restarts.jl b/test/unit_tests/test_restarts.jl index 5b060ae2..4c754d4a 100644 --- a/test/unit_tests/test_restarts.jl +++ b/test/unit_tests/test_restarts.jl @@ -2,7 +2,6 @@ test_path = joinpath(HallThruster.PACKAGE_ROOT, "test", "unit_tests") restart_path = joinpath(test_path, "test_restart.jld2") - config = HallThruster.Config(; ncharge = 1, anode_Te = 2.0, @@ -28,12 +27,10 @@ grid = HallThruster.generate_grid(config.thruster.geometry, config.domain, EvenGrid(ncells)) fluids, fluid_ranges, species, species_range_dict = HallThruster.configure_fluids(config) - - #make a new restart from the current configuration - sol = HallThruster.run_simulation(config; dt = 5e-9, duration=4e-9, grid = HallThruster.EvenGrid(ncells), nsave = 10) + #make a new restart from the current configuration + sol = HallThruster.run_simulation(config; dt = 1e-9, duration=4e-9, grid = HallThruster.EvenGrid(ncells), nsave = 10) HallThruster.write_restart(restart_path, sol) - # Check that writing a restart of a restart does not lose any information # i.e. restarts are perfectly reconstructable restart = HallThruster.read_restart(restart_path) @@ -59,8 +56,8 @@ @test cache.Tev ≈ sol[:Tev][end] @test cache.νc ≈ sol[:νc][end] - - #test that a simulation can be run from the restart + + #test that a simulation can be run from the restart restart_sol = HallThruster.run_simulation(config; dt = 5e-9, duration=4e-9, grid = HallThruster.EvenGrid(ncells), nsave = 10, restart = restart_path) @test restart_sol.retcode == :success From 1b3cc5c58244df5cbbce16915b3ab82b76251ec7 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 22:30:04 -0400 Subject: [PATCH 07/14] tests pass --- src/HallThruster.jl | 1 + src/simulation/solution.jl | 2 +- src/simulation/update_heavy_species.jl | 6 +++--- test/order_verification/ovs_energy.jl | 5 +++-- test/order_verification/ovs_ions.jl | 5 ++--- test/runtests.jl | 1 + test/unit_tests/test_restarts.jl | 28 ++++++++------------------ 7 files changed, 19 insertions(+), 29 deletions(-) diff --git a/src/HallThruster.jl b/src/HallThruster.jl index 268962c6..846e4eb5 100644 --- a/src/HallThruster.jl +++ b/src/HallThruster.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 diff --git a/src/simulation/solution.jl b/src/simulation/solution.jl index 46997605..03e8170b 100644 --- a/src/simulation/solution.jl +++ b/src/simulation/solution.jl @@ -54,7 +54,7 @@ function solve(U, params, tspan; saveat) integrate_heavy_species!(U, params, params.dt[]) update_heavy_species!(U, params) - # Check for NaNs and terminate if necessary + # Check for NaNs in heavy species solve and terminate if necessary nandetected = false infdetected = false diff --git a/src/simulation/update_heavy_species.jl b/src/simulation/update_heavy_species.jl index 05b94c31..6d37b413 100644 --- a/src/simulation/update_heavy_species.jl +++ b/src/simulation/update_heavy_species.jl @@ -66,16 +66,16 @@ function iterate_heavy_species!(dU, U, params; apply_boundary_conditions = true) end # Perform one step of the Strong-stability-preserving RK22 algorithm to the ion fluid -function integrate_heavy_species!(U, params, dt) +function integrate_heavy_species!(U, params, dt, apply_boundary_conditions = true) (;k, u1) = params.cache # First step of SSPRK22 - iterate_heavy_species!(k, U, params) + iterate_heavy_species!(k, U, params; apply_boundary_conditions) @. u1 = U + dt * k stage_limiter!(u1, params) # Second step of SSPRK22 - iterate_heavy_species!(k, u1, params) + iterate_heavy_species!(k, u1, params; apply_boundary_conditions) @. U = (U + u1 + dt * k) / 2 stage_limiter!(U, params) diff --git a/test/order_verification/ovs_energy.jl b/test/order_verification/ovs_energy.jl index 42ad5c91..af3da07a 100644 --- a/test/order_verification/ovs_energy.jl +++ b/test/order_verification/ovs_energy.jl @@ -151,7 +151,8 @@ function verify_energy(ncells; niters = 20000) ionization_product_indices, excitation_reactions, excitation_reactant_indices, - Δz_cell, Δz_edge + Δz_cell, Δz_edge, + ncells ) solve_energy!(params, niters, dt) @@ -177,7 +178,7 @@ function verify_energy(ncells; niters = 20000) ionization_product_indices, excitation_reactions, excitation_reactant_indices, - Δz_cell, Δz_edge + Δz_cell, Δz_edge, ncells ) solve_energy!(params, niters, dt) diff --git a/test/order_verification/ovs_ions.jl b/test/order_verification/ovs_ions.jl index b008962e..a9b82152 100644 --- a/test/order_verification/ovs_ions.jl +++ b/test/order_verification/ovs_ions.jl @@ -114,8 +114,6 @@ function solve_ions(ncells, scheme; t_end = 1e-4) ionization_reactions = HallThruster._load_reactions(config.ionization_model, species) index = HallThruster.configure_index(fluids, fluid_ranges) - @show is_velocity_index, fluid_ranges - nvars = fluid_ranges[end][end] F = zeros(nvars, nedges) @@ -150,6 +148,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) cache = (;k, u1, ue, μ, F, UL, UR, ∇ϕ, λ_global, channel_area, dA_dz, dt_cell, dt, dt_u, dt_iz, dt_E, nϵ, νiz, inelastic_losses) params = (; + ncells = ncells+2, index, config, cache, @@ -176,7 +175,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) t = 0.0 while t < t_end @views U[:, end] = U[:, end-1] - HallThruster.integrate_heavy_species!(U, params, dt[]); + HallThruster.integrate_heavy_species!(U, params, dt[], false); t += dt[] end diff --git a/test/runtests.jl b/test/runtests.jl index 9e377701..b867ab52 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ using SparseArrays doctest(HallThruster) +# exercise all parts of the solver loop HallThruster.example_simulation(;ncells=20, duration=1e-7, dt=1e-8, nsave=2) include("unit_tests/test_gas.jl") diff --git a/test/unit_tests/test_restarts.jl b/test/unit_tests/test_restarts.jl index 4c754d4a..4b06814a 100644 --- a/test/unit_tests/test_restarts.jl +++ b/test/unit_tests/test_restarts.jl @@ -3,24 +3,11 @@ restart_path = joinpath(test_path, "test_restart.jld2") config = HallThruster.Config(; - ncharge = 1, - anode_Te = 2.0, - cathode_Te = 2.0, - discharge_voltage = 300.0, - excitation_model = HallThruster.LandmarkExcitationLookup(), - wall_loss_model = HallThruster.WallSheath(HallThruster.BoronNitride, 1.0), - implicit_energy = 1.0, - neutral_velocity = 300.0, - neutral_temperature = 300.0, - ion_temperature = 1000.0, thruster = HallThruster.SPT_100, - anode_mass_flow_rate = 5e-6, - electron_neutral_model = HallThruster.LandmarkElectronNeutral(), - electron_ion_collisions = true, - ionization_model = HallThruster.LandmarkIonizationLookup(), - domain = (0.0, 0.05), - LANDMARK = true, - conductivity_model = HallThruster.LANDMARK_conductivity(), + domain = (0.0u"cm", 8.0u"cm"), + discharge_voltage = 300.0u"V", + anode_mass_flow_rate = 5u"mg/s", + wall_loss_model = HallThruster.WallSheath(HallThruster.BoronNitride, 0.15) ) ncells = 50 @@ -28,7 +15,10 @@ fluids, fluid_ranges, species, species_range_dict = HallThruster.configure_fluids(config) #make a new restart from the current configuration - sol = HallThruster.run_simulation(config; dt = 1e-9, duration=4e-9, grid = HallThruster.EvenGrid(ncells), nsave = 10) + duration=1e-7 + dt=1e-8 + nsave=2 + sol = HallThruster.run_simulation(config; ncells, duration, dt, nsave, verbose = false) HallThruster.write_restart(restart_path, sol) # Check that writing a restart of a restart does not lose any information @@ -56,7 +46,6 @@ @test cache.Tev ≈ sol[:Tev][end] @test cache.νc ≈ sol[:νc][end] - #test that a simulation can be run from the restart restart_sol = HallThruster.run_simulation(config; dt = 5e-9, duration=4e-9, grid = HallThruster.EvenGrid(ncells), nsave = 10, restart = restart_path) @test restart_sol.retcode == :success @@ -95,5 +84,4 @@ end @test landmark_1[3].u[1][1, :] ≈ landmark_1_hybrid.u[1][1, :] @test landmark_1[3].u[1][2, :] ≈ landmark_1_hybrid.u[1][2, :] @test landmark_1[3].u[1][4, :] ≈ landmark_1_hybrid.u[1][4, :] - end From 6c8bc718b1f9ebca718983d8c8fd75ab4882a51d Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 22:34:26 -0400 Subject: [PATCH 08/14] remove solve_background_neutrals flag --- docs/src/config.md | 3 +-- src/simulation/boundaryconditions.jl | 4 +--- src/simulation/configuration.jl | 3 --- src/simulation/simulation.jl | 12 +++--------- test/order_verification/ovs_energy.jl | 4 ++-- test/order_verification/ovs_ions.jl | 1 - test/unit_tests/input_shifted.json | 1 - test/unit_tests/input_twozone.json | 1 - test/unit_tests/test_boundary_conditions.jl | 1 - test/unit_tests/test_initialization.jl | 2 -- test/unit_tests/test_json.jl | 3 +-- test/unit_tests/test_misc.jl | 11 ----------- 12 files changed, 8 insertions(+), 38 deletions(-) diff --git a/docs/src/config.md b/docs/src/config.md index 60b9c452..3e7f8254 100644 --- a/docs/src/config.md +++ b/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 diff --git a/src/simulation/boundaryconditions.jl b/src/simulation/boundaryconditions.jl index d35bf608..f9a68b38 100644 --- a/src/simulation/boundaryconditions.jl +++ b/src/simulation/boundaryconditions.jl @@ -29,9 +29,7 @@ function left_boundary_state!(bc_state, U, params) 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] += params.background_neutral_density * params.background_neutral_velocity / un - end + bc_state[index.ρn] += params.background_neutral_density * params.background_neutral_velocity / un @inbounds for Z in 1:params.config.ncharge interior_density = U[index.ρi[Z], begin+1] diff --git a/src/simulation/configuration.jl b/src/simulation/configuration.jl index 0ae1f964..12f7792e 100644 --- a/src/simulation/configuration.jl +++ b/src/simulation/configuration.jl @@ -41,7 +41,6 @@ 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 anode_boundary_condition::Symbol @@ -85,7 +84,6 @@ 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", anode_boundary_condition = :sheath, @@ -146,7 +144,6 @@ function Config(; LANDMARK, anode_mass_flow_rate, ion_wall_losses, - solve_background_neutrals, background_pressure, background_neutral_temperature, anode_boundary_condition, diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 96589c3b..7b65a1fa 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -132,11 +132,6 @@ function run_simulation( error("LANDMARK configuration needs to use the LANDMARK thermal conductivity model.") end - #check that user is aware of neutral backflow flag - if (config.background_pressure > 0.0) && (config.solve_background_neutrals == false) - @warn("Background neutral pressure set but solve background neutrals not enabled. Did you mean to set solve_background_neutrals to true?") - end - # If duration and/or dt are provided with units, convert to seconds and then strip units duration = convert_to_float64(duration, u"s") dt = convert_to_float64(dt, u"s") @@ -291,7 +286,7 @@ function run_simulation(json_content::JSON3.Object; single_section = false, nons anom_model, cathode_location_m, max_charge, num_cells, dt_s, duration_s, num_save, flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, solve_background_neutrals, + ion_wall_losses, electron_ion_collisions, # Parameters sheath_loss_coefficient, ion_temp_K, neutral_temp_K, neutral_velocity_m_s, @@ -313,7 +308,7 @@ function run_simulation(json_content::JSON3.Object; single_section = false, nons anom_model, cathode_location_m, max_charge, num_cells, dt_s, duration_s, num_save, flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, solve_background_neutrals, + ion_wall_losses, electron_ion_collisions, # Parameters sheath_loss_coefficient, ion_temp_K, neutral_temp_K, neutral_velocity_m_s, @@ -358,7 +353,7 @@ function run_simulation(json_content::JSON3.Object; single_section = false, nons anom_model, cathode_location_m, max_charge, num_cells, dt_s, duration_s, num_save, flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, solve_background_neutrals, + ion_wall_losses, electron_ion_collisions, ) = simulation (; @@ -444,7 +439,6 @@ function run_simulation(json_content::JSON3.Object; single_section = false, nons scheme = HyperbolicScheme(; flux_function = eval(Symbol(flux_function)), limiter = eval(Symbol(limiter)), reconstruct ), - solve_background_neutrals = solve_background_neutrals, background_pressure = background_pressure_Torr * u"Torr", background_neutral_temperature = background_temperature_K * u"K", apply_thrust_divergence_correction diff --git a/test/order_verification/ovs_energy.jl b/test/order_verification/ovs_energy.jl index af3da07a..b148d1c1 100644 --- a/test/order_verification/ovs_energy.jl +++ b/test/order_verification/ovs_energy.jl @@ -118,7 +118,7 @@ function verify_energy(ncells; niters = 20000) ncharge = 1, source_energy = source_func, implicit_energy = 1.0, min_electron_temperature, transition_function, LANDMARK, propellant, ionization_model, excitation_model, wall_loss_model, geometry, - anode_boundary_condition = :dirichlet, solve_background_neutrals = false, + anode_boundary_condition = :dirichlet, conductivity_model = HallThruster.LANDMARK_conductivity(), ) @@ -165,7 +165,7 @@ function verify_energy(ncells; niters = 20000) ncharge = 1, source_energy = source_func, implicit_energy = 0.5, min_electron_temperature, transition_function, LANDMARK, propellant, ionization_model, excitation_model, wall_loss_model, geometry, - anode_boundary_condition = :dirichlet, solve_background_neutrals = false, + anode_boundary_condition = :dirichlet, conductivity_model = HallThruster.LANDMARK_conductivity(), ) diff --git a/test/order_verification/ovs_ions.jl b/test/order_verification/ovs_ions.jl index a9b82152..ab5cecb3 100644 --- a/test/order_verification/ovs_ions.jl +++ b/test/order_verification/ovs_ions.jl @@ -94,7 +94,6 @@ function solve_ions(ncells, scheme; t_end = 1e-4) conductivity_model = HallThruster.LANDMARK_conductivity(), ion_wall_losses = false, anode_boundary_condition = :dirichlet, - solve_background_neutrals = false, ) z_edge = grid.edges diff --git a/test/unit_tests/input_shifted.json b/test/unit_tests/input_shifted.json index 3fdf2ee7..21da469a 100644 --- a/test/unit_tests/input_shifted.json +++ b/test/unit_tests/input_shifted.json @@ -42,7 +42,6 @@ "ion_wall_losses": true, "electron_ion_collisions": true, "anom_model": "ShiftedTwoZone", - "solve_background_neutrals": true, "adaptive": true } } \ No newline at end of file diff --git a/test/unit_tests/input_twozone.json b/test/unit_tests/input_twozone.json index 80599ae5..9240402f 100644 --- a/test/unit_tests/input_twozone.json +++ b/test/unit_tests/input_twozone.json @@ -38,7 +38,6 @@ "ion_wall_losses": true, "electron_ion_collisions": true, "anom_model": "TwoZoneBohm", - "solve_background_neutrals": true, "adaptive": true } } \ No newline at end of file diff --git a/test/unit_tests/test_boundary_conditions.jl b/test/unit_tests/test_boundary_conditions.jl index 5c8588e3..d80b6675 100644 --- a/test/unit_tests/test_boundary_conditions.jl +++ b/test/unit_tests/test_boundary_conditions.jl @@ -13,7 +13,6 @@ LANDMARK = true, conductivity_model = HallThruster.LANDMARK_conductivity(), anode_boundary_condition = :dirichlet, - solve_background_neutrals = :true, neutral_temperature = 300.0, ion_temperature = 300.0, background_pressure = 6e-4, diff --git a/test/unit_tests/test_initialization.jl b/test/unit_tests/test_initialization.jl index 337341c5..40a93b7c 100644 --- a/test/unit_tests/test_initialization.jl +++ b/test/unit_tests/test_initialization.jl @@ -103,7 +103,6 @@ end ) config = HallThruster.Config(; - solve_background_neutrals = false, background_pressure = 0.0u"Torr", background_neutral_temperature = 0.0u"K", common_opts... @@ -149,7 +148,6 @@ end TB = 120u"K" config_bg = HallThruster.Config(; - solve_background_neutrals = true, background_pressure = pB, background_neutral_temperature = TB, common_opts... diff --git a/test/unit_tests/test_json.jl b/test/unit_tests/test_json.jl index d2fd3d59..3292675f 100644 --- a/test/unit_tests/test_json.jl +++ b/test/unit_tests/test_json.jl @@ -15,7 +15,6 @@ @test config.thruster.name == "SPT-100" @test config.propellant == Xenon @test config.anode_mass_flow_rate ≈ 3e-6 - @test config.solve_background_neutrals == true @test config.ion_wall_losses == true @test sol.params.adaptive == true @@ -25,4 +24,4 @@ @test config.anom_model isa HallThruster.TwoZoneBohm @test config.anom_model.coeffs[1] ≈ 1/160 @test config.anom_model.coeffs[2] ≈ 1/16 -end \ No newline at end of file +end diff --git a/test/unit_tests/test_misc.jl b/test/unit_tests/test_misc.jl index 979a58c8..d75e2760 100644 --- a/test/unit_tests/test_misc.jl +++ b/test/unit_tests/test_misc.jl @@ -17,17 +17,6 @@ ) @test_logs (:warn, "CFL for adaptive timestepping set higher than stability limit of 0.8. Setting CFL to 0.799.") HallThruster.run_simulation(config; dt=5e-9, duration=0e-9, grid = HallThruster.EvenGrid(2), nsave = 10, adaptive = true, CFL = 0.9) - - pressure_config = HallThruster.Config(; - thruster = HallThruster.SPT_100, - domain = (0.0u"cm", 8.0u"cm"), - discharge_voltage = 300.0u"V", - anode_mass_flow_rate = 5u"mg/s", - background_pressure = 1.0u"Pa", - ) - - @test_logs (:warn, "Background neutral pressure set but solve background neutrals not enabled. Did you mean to set solve_background_neutrals to true?") HallThruster.run_simulation(pressure_config; dt=5e-9, duration=4e-9, grid = HallThruster.EvenGrid(2), nsave = 10) - end @testset "Linear algebra" begin From 770fc24efe09f4d7bd2c949688fbd41b512c2430 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 24 May 2024 23:16:19 -0400 Subject: [PATCH 09/14] update version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 85bc1310..d8d516c0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "HallThruster" uuid = "2311f341-5e6d-4941-9e3e-3ce0ae0d9ed6" authors = ["Thomas Marks "] -version = "0.13.2" +version = "0.14.0" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" From 3faf4cb3eb2c699238a5838937f332a202046ec9 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Wed, 5 Jun 2024 13:31:19 -0400 Subject: [PATCH 10/14] update OVS tests --- test/order_verification/ovs_ions.jl | 17 ++++------------- test/runtests.jl | 11 ++++------- 2 files changed, 8 insertions(+), 20 deletions(-) diff --git a/test/order_verification/ovs_ions.jl b/test/order_verification/ovs_ions.jl index ab5cecb3..90740dbd 100644 --- a/test/order_verification/ovs_ions.jl +++ b/test/order_verification/ovs_ions.jl @@ -18,29 +18,22 @@ const e = HallThruster.e const Ti = 300 const L = 0.05 -ϕ = 0.0 * sin_wave(x/L, amplitude = 300, phase = π/2, nwaves = 0.5) +ϕ = sin_wave(x/L, amplitude = 5, phase = π/2, nwaves = 0.5) ne = sin_wave(x/L, amplitude = 1e13, phase = π/2, nwaves = 2, offset = 6e13) nn = sin_wave(x/L, amplitude = 2e18, phase = π/2, nwaves = 1, offset = 6e18) ui = sin_wave(x/L, amplitude = 2000, phase = -π/2, nwaves = 0.5, offset = 3000) -μ = sin_wave(x/L, amplitude = 1e2, phase = 3π/2, nwaves = 0.6, offset = 1.1e2) ϵ = sin_wave(x/L, amplitude = 3, phase = -π/2, nwaves = 1, offset = 6) nϵ = ne * ϵ -pe = e * nϵ ∇ϕ = Dx(ϕ) -∇pe = Dx(pe) ρi = ne * mi ρiui = ρi * ui ρn = nn * mi -ue = μ * (∇ϕ - Dx(nϵ)/ne) p = ne * HallThruster.kB * Ti -∇p = Dx(p) ϕ_func = eval(build_function(ϕ, [x])) ne_func = eval(build_function(ne, [x])) -μ_func = eval(build_function(μ, [x])) ρiui_func = eval(build_function(ρiui, [x])) nϵ_func = eval(build_function(nϵ, [x])) -ue_func = eval(build_function(expand_derivatives(ue), [x])) ∇ϕ_func = eval(build_function(expand_derivatives(∇ϕ), [x])) ρn_func = eval(build_function(ρn, [x])) ρi_func = eval(build_function(ρi, [x])) @@ -101,8 +94,6 @@ function solve_ions(ncells, scheme; t_end = 1e-4) nedges = length(z_edge) - ue = ue_func.(z_cell) - μ = μ_func.(z_cell) nϵ = nϵ_func.(z_cell) ρn_exact = ρn_func.(z_cell) ρi_exact = ρi_func.(z_cell) @@ -144,7 +135,7 @@ function solve_ions(ncells, scheme; t_end = 1e-4) u1 = zeros(size(U)) inelastic_losses = zeros(ncells+2) νiz = zeros(ncells+2) - cache = (;k, u1, ue, μ, F, UL, UR, ∇ϕ, λ_global, channel_area, dA_dz, dt_cell, dt, dt_u, dt_iz, dt_E, nϵ, νiz, inelastic_losses) + cache = (;k, u1, F, UL, UR, ∇ϕ, λ_global, channel_area, dA_dz, dt_cell, dt, dt_u, dt_iz, dt_E, nϵ, νiz, inelastic_losses) params = (; ncells = ncells+2, @@ -174,8 +165,8 @@ function solve_ions(ncells, scheme; t_end = 1e-4) t = 0.0 while t < t_end @views U[:, end] = U[:, end-1] - HallThruster.integrate_heavy_species!(U, params, dt[], false); - t += dt[] + HallThruster.integrate_heavy_species!(U, params, cache.dt[], false); + t += cache.dt[] end sol = (;t = [t], u = [U]) diff --git a/test/runtests.jl b/test/runtests.jl index b867ab52..a8777e92 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -51,12 +51,13 @@ end include("order_verification/ovs_ions.jl") @testset "Order verification (neutrals and ions)" begin - refinements = refines(6, 10, 2) + refinements = refines(5, 10, 2) limiter = HallThruster.van_leer flux_names = ["HLLE", "Rusanov", "Global Lax-Friedrichs"] #flux_names = ["Global Lax-Friedrichs"] fluxes = [HallThruster.HLLE, HallThruster.rusanov, HallThruster.global_lax_friedrichs,] + #fluxes = [HallThruster.HLLE, HallThruster.rusanov, HallThruster.global_lax_friedrichs,] # Which L-P norms to check cases = ["ρn", "ρi", "ρiui"] @@ -65,13 +66,9 @@ include("order_verification/ovs_ions.jl") for (flux, flux_name) in zip(fluxes, flux_names) for reconstruct in [false, true] - if reconstruct && flux_name == "HLLE" - continue - end scheme = HallThruster.HyperbolicScheme(flux, limiter, reconstruct) orders, norms = test_refinements(ncells -> OVS_Ions.solve_ions(ncells, scheme), refinements, norms_to_test) for (i, (order, norm)) in enumerate(zip(orders, norms)) - norm_ind = mod1(i, num_norms) case_ind = ((i-1) ÷ num_norms) + 1 case_str = "($(cases[case_ind]), $(norms_to_test[norm_ind])-norm)" @@ -84,9 +81,9 @@ include("order_verification/ovs_ions.jl") # Check that we achieve the desired order of accuracy if (reconstruct) - @test order > 1.5 + @test order >= 1.5 else - @test order > 0.85 + @test order > 0.75 end end end From 4d482051bc5aa54f5479510acf4204efc722d47a Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Wed, 5 Jun 2024 13:50:26 -0400 Subject: [PATCH 11/14] remove transition functions and update wall losses --- src/HallThruster.jl | 1 - src/collisions/anomalous.jl | 4 +- src/simulation/configuration.jl | 15 +++-- src/simulation/initialization.jl | 2 +- src/simulation/simulation.jl | 2 +- src/simulation/sourceterms.jl | 5 +- src/simulation/update_electrons.jl | 3 +- src/utilities/transition_functions.jl | 65 ------------------- src/utilities/utility_functions.jl | 14 ++++ .../constant_sheath_potential.jl | 2 +- src/wall_loss_models/wall_sheath.jl | 28 +++++--- test/order_verification/ovs_energy.jl | 6 +- test/unit_tests/test_electrons.jl | 10 +-- test/unit_tests/test_misc.jl | 33 ---------- test/unit_tests/test_walls.jl | 25 +++---- 15 files changed, 73 insertions(+), 142 deletions(-) delete mode 100644 src/utilities/transition_functions.jl diff --git a/src/HallThruster.jl b/src/HallThruster.jl index 846e4eb5..675e5995 100644 --- a/src/HallThruster.jl +++ b/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") diff --git a/src/collisions/anomalous.jl b/src/collisions/anomalous.jl index 5c00c3cf..e671d281 100644 --- a/src/collisions/anomalous.jl +++ b/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 diff --git a/src/simulation/configuration.jl b/src/simulation/configuration.jl index 12f7792e..2f088e8e 100644 --- a/src/simulation/configuration.jl +++ b/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 @@ -43,6 +43,7 @@ struct Config{A<:AnomalousTransportModel, TC<:ThermalConductivityModel, W<:WallL ion_wall_losses::Bool background_pressure::Float64 background_neutral_temperature::Float64 + neutral_ingestion_multiplier::Float64 anode_boundary_condition::Symbol anom_smoothing_iters::Int solve_plume::Bool @@ -58,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, @@ -73,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, @@ -86,6 +87,7 @@ function Config(; ion_wall_losses = 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, @@ -123,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 @@ -132,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, @@ -146,6 +150,7 @@ function Config(; ion_wall_losses, background_pressure, background_neutral_temperature, + neutral_ingestion_multiplier, anode_boundary_condition, anom_smoothing_iters, solve_plume, diff --git a/src/simulation/initialization.jl b/src/simulation/initialization.jl index 7ea37476..3b560d77 100644 --- a/src/simulation/initialization.jl +++ b/src/simulation/initialization.jl @@ -44,7 +44,7 @@ function initialize!(U, params, ::DefaultInitialization) # Beam neutral density at outlet ρn_1 = 0.01 * ρn_0 - neutral_function = z -> SmoothIf(transition_length = L_ch / 6)(z-z0, L_ch / 2, ρn_0, ρn_1) + neutral_function = z -> smooth_if(z-z0, L_ch / 2, ρn_0, ρn_1, L_ch / 6) number_density_function = z -> sum(Z * ion_density_function(z, Z) / mi for Z in 1:ncharge) diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 7b65a1fa..92262497 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -435,7 +435,7 @@ function run_simulation(json_content::JSON3.Object; single_section = false, nons neutral_velocity = neutral_velocity_m_s, electron_ion_collisions = electron_ion_collisions, min_electron_temperature = cathode_electron_temp_eV, - transition_function = LinearTransition(inner_outer_transition_length_m, 0.0), + transition_length = inner_outer_transition_length_m, scheme = HyperbolicScheme(; flux_function = eval(Symbol(flux_function)), limiter = eval(Symbol(limiter)), reconstruct ), diff --git a/src/simulation/sourceterms.jl b/src/simulation/sourceterms.jl index 5e212234..32d9e49c 100644 --- a/src/simulation/sourceterms.jl +++ b/src/simulation/sourceterms.jl @@ -85,9 +85,10 @@ function apply_ion_wall_losses!(dU, U, params, i) end @inbounds for Z in 1:ncharge - in_channel = config.transition_function(z_cell[i], L_ch, 1.0, 0.0) + in_channel = linear_transition(z_cell[i], L_ch, params.config.transition_length, 1.0, 0.0) u_bohm = sqrt(Z * e * params.cache.Tev[i] / mi) - νiw = α * in_channel * u_bohm / Δr + h = edge_to_center_density_ratio() + νiw = α * in_channel * u_bohm / Δr * h density_loss = U[index.ρi[Z], i] * νiw momentum_loss = U[index.ρiui[Z], i] * νiw diff --git a/src/simulation/update_electrons.jl b/src/simulation/update_electrons.jl index ffa0a244..220e40de 100644 --- a/src/simulation/update_electrons.jl +++ b/src/simulation/update_electrons.jl @@ -35,8 +35,7 @@ function update_electrons!(params, t = 0) # Compute wall collision frequency, with transition function to force no momentum wall collisions in plume radial_loss_frequency[i] = freq_electron_wall(params.config.wall_loss_model, params, i) - νew_momentum[i] = radial_loss_frequency[i]* params.config.transition_function(params.z_cell[i], params.L_ch, 1.0, 0.0) - + νew_momentum[i] = radial_loss_frequency[i]* linear_transition(params.z_cell[i], params.L_ch, params.config.transition_length, 1.0, 0.0) end # Update anomalous transport diff --git a/src/utilities/transition_functions.jl b/src/utilities/transition_functions.jl deleted file mode 100644 index 8f2d371f..00000000 --- a/src/utilities/transition_functions.jl +++ /dev/null @@ -1,65 +0,0 @@ -abstract type TransitionFunction end - -Base.@kwdef struct SmoothIf <: TransitionFunction - transition_length::Float64 -end - -@inline smooth_if(x, cutoff, y1, y2, L) = ((y2 - y1)*tanh((x-cutoff)/(L/4)) + y1 + y2) / 2 - -(s::SmoothIf)(x, cutoff, y1, y2) = smooth_if(x, cutoff, y1, y2, s.transition_length) - -Base.@kwdef struct SmoothStep <: TransitionFunction - transition_length::Float64 -end - -@inline smoothstep(x, x1, x2, y1, y2) = let t = (x - x1) / (x2 - x1) - y = if t < 0 - 0.0 - elseif t > 1 - 1.0 - else - 3 * t^2 - 2 * t^3 - end - return (y2 - y1) * y + y1 -end - -(s::SmoothStep)(x, cutoff, y1, y2) = let L = s.transition_length - smoothstep(x, cutoff - L/2, cutoff + L/2, y1, y2) -end - -Base.@kwdef struct QuadraticTransition <: TransitionFunction - transition_length::Float64 - offset::Float64 = 0.0 -end - -function (q::QuadraticTransition)(x, cutoff, y1, y2) - x′ = x - q.offset - if x′ < cutoff - return y1 - else - return y1 + (y2 - y1) * (x′ - cutoff)^2 / ((x - cutoff)^2 + q.transition_length^2) - end -end - -Base.@kwdef struct LinearTransition <: TransitionFunction - transition_length::Float64 - offset::Float64 = 0.0 -end - -function(ℓ::LinearTransition)(x, cutoff, y1, y2) - x′ = x - ℓ.offset - L = ℓ.transition_length - x1 = cutoff - L/2 - x2 = cutoff + L/2 - if x′ < x1 - return y1 - elseif x′ > x2 - return y2 - else - return lerp(x, x1, x2, y1, y2) - end -end - -struct StepFunction <: TransitionFunction end - -(::StepFunction)(x, cutoff, y1, y2) = x ≤ cutoff ? y1 : y2 diff --git a/src/utilities/utility_functions.jl b/src/utilities/utility_functions.jl index d9c071f9..6b7990c4 100644 --- a/src/utilities/utility_functions.jl +++ b/src/utilities/utility_functions.jl @@ -56,6 +56,20 @@ julia> lerp(0.5, 0.0, 1.0, 0.0, 2.0) return muladd(t, (y1 - y0), y0) end +@inline function linear_transition(x, cutoff, L, y1, y2) + x1 = cutoff - L/2 + x2 = cutoff + L/2 + if x < x1 + return y1 + elseif x > x2 + return y2 + else + return lerp(x, x1, x2, y1, y2) + end +end + +@inline smooth_if(x, cutoff, y1, y2, L) = ((y2 - y1)*tanh((x-cutoff)/(L/4)) + y1 + y2) / 2 + function find_left_index(value, array) N = length(array) diff --git a/src/wall_loss_models/constant_sheath_potential.jl b/src/wall_loss_models/constant_sheath_potential.jl index 7f226aa0..b3fdee8b 100644 --- a/src/wall_loss_models/constant_sheath_potential.jl +++ b/src/wall_loss_models/constant_sheath_potential.jl @@ -14,7 +14,7 @@ function wall_power_loss(model::ConstantSheathPotential, params, i) ϵ = nϵ / ne (;sheath_potential, inner_loss_coeff, outer_loss_coeff) = model - αϵ = config.transition_function(z_cell[i], L_ch, inner_loss_coeff, outer_loss_coeff) + αϵ = linear_transition(z_cell[i], L_ch, config.transition_length, inner_loss_coeff, outer_loss_coeff) W = 1e7 * αϵ * ϵ * exp(-sheath_potential / ϵ) diff --git a/src/wall_loss_models/wall_sheath.jl b/src/wall_loss_models/wall_sheath.jl index 624c8153..201450bb 100644 --- a/src/wall_loss_models/wall_sheath.jl +++ b/src/wall_loss_models/wall_sheath.jl @@ -28,7 +28,7 @@ function wall_electron_temperature(params, i) L_ch = config.thruster.geometry.channel_length - Tev = config.transition_function(z_cell[i], L_ch, Tev_channel, Tev_plume) + Tev = linear_transition(z_cell[i], L_ch, config.transition_length, Tev_channel, Tev_plume) return Tev end @@ -49,8 +49,13 @@ Base.@kwdef struct WallSheath <: WallLossModel end end +# compute the edge-to-center density ratio +# (c.f https://iopscience.iop.org/article/10.1088/0963-0252/24/2/025017) +# this is approximate, but deviations only become noticable when the knudsen number becomes small, which is not true in our case +@inline edge_to_center_density_ratio() = 0.86 / sqrt(3); + function freq_electron_wall(model::WallSheath, params, i) - (; index, config, cache) = params + (; config, cache) = params (; ncharge) = config mi = config.propellant.m #compute radii difference @@ -61,20 +66,23 @@ function freq_electron_wall(model::WallSheath, params, i) #calculate and store SEE coefficient γ = SEE_yield(model.material, Tev, params.γ_SEE_max) cache.γ_SEE[i] = γ - #compute the ion current to the walls + + # compute the ion current to the walls + h = edge_to_center_density_ratio() j_iw = 0.0 for Z in 1:ncharge - niw = cache.ni[Z, i] - j_iw += model.α * Z * niw * sqrt(Z * e * Tev / mi) + niw = h * cache.ni[Z, i] + j_iw += Z * model.α * niw * sqrt(Z * e * Tev / mi) end - #compute electron wall collision frequency + + # compute electron wall collision frequency νew = j_iw / (Δr * (1 - γ)) / cache.ne[i] return νew end function wall_power_loss(model::WallSheath, params, i) - (;config) = params + (;config, cache, z_cell, L_ch) = params mi = config.propellant.m Tev = wall_electron_temperature(params, i) @@ -86,10 +94,12 @@ function wall_power_loss(model::WallSheath, params, i) ϕ_s = sheath_potential(Tev, γ, mi) # Compute electron wall collision frequency with transition function for energy wall collisions in plume - νew = params.cache.radial_loss_frequency[i] * params.config.transition_function(params.z_cell[i], params.L_ch, 1.0, params.config.electron_plume_loss_scale) + νew = cache.radial_loss_frequency[i] * linear_transition( + z_cell[i], L_ch, config.transition_length, 1.0, config.electron_plume_loss_scale + ) # Compute wall power loss rate - W = νew * (2 * (1 - 0.5 * model.material.σ₀) * Tev + (1 - γ) * ϕ_s) / γ + W = νew * (2 * Tev + (1 - γ) * ϕ_s) return W end diff --git a/test/order_verification/ovs_energy.jl b/test/order_verification/ovs_energy.jl index b148d1c1..1fb1d550 100644 --- a/test/order_verification/ovs_energy.jl +++ b/test/order_verification/ovs_energy.jl @@ -102,7 +102,7 @@ function verify_energy(ncells; niters = 20000) # Test backward difference implicit solve dt = 1e-6 - transition_function = HallThruster.StepFunction() + transition_length = 0.0 excitation_model = OVS_Excitation() ionization_model = OVS_Ionization() @@ -116,7 +116,7 @@ function verify_energy(ncells; niters = 20000) config = (; ncharge = 1, source_energy = source_func, implicit_energy = 1.0, - min_electron_temperature, transition_function, LANDMARK, propellant, + min_electron_temperature, transition_length, LANDMARK, propellant, ionization_model, excitation_model, wall_loss_model, geometry, anode_boundary_condition = :dirichlet, conductivity_model = HallThruster.LANDMARK_conductivity(), @@ -163,7 +163,7 @@ function verify_energy(ncells; niters = 20000) config = (; ncharge = 1, source_energy = source_func, implicit_energy = 0.5, - min_electron_temperature, transition_function, LANDMARK, propellant, + min_electron_temperature, transition_length, LANDMARK, propellant, ionization_model, excitation_model, wall_loss_model, geometry, anode_boundary_condition = :dirichlet, conductivity_model = HallThruster.LANDMARK_conductivity(), diff --git a/test/unit_tests/test_electrons.jl b/test/unit_tests/test_electrons.jl index d57dda4c..6cb77d8f 100644 --- a/test/unit_tests/test_electrons.jl +++ b/test/unit_tests/test_electrons.jl @@ -29,14 +29,14 @@ c2 = 1/16 anom_model = HallThruster.TwoZoneBohm(c1, c2) thruster = HallThruster.SPT_100 - transition_function = HallThruster.StepFunction() + transition_length = 0.0 wall_collision_freq = 1e7 - config_landmark = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.LandmarkElectronNeutral(), electron_ion_collisions = false, ncharge = 1, thruster, transition_function, wall_collision_freq = 1e7) - config_gk = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.GKElectronNeutral(), electron_ion_collisions = true, ncharge = 1, thruster, transition_function, wall_collision_freq = 1e7) - config_complex = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.ElectronNeutralLookup(), electron_ion_collisions = true, ncharge = 1, thruster, transition_function, wall_collision_freq = 1e7) - config_none = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.NoElectronNeutral(), electron_ion_collisions = false, ncharge = 1, thruster, transition_function, wall_collision_freq = 1e7) + config_landmark = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.LandmarkElectronNeutral(), electron_ion_collisions = false, ncharge = 1, thruster, transition_length, wall_collision_freq = 1e7) + config_gk = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.GKElectronNeutral(), electron_ion_collisions = true, ncharge = 1, thruster, transition_length, wall_collision_freq = 1e7) + config_complex = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.ElectronNeutralLookup(), electron_ion_collisions = true, ncharge = 1, thruster, transition_length, wall_collision_freq = 1e7) + config_none = (;anom_model, propellant = HallThruster.Xenon, electron_neutral_model = HallThruster.NoElectronNeutral(), electron_ion_collisions = false, ncharge = 1, thruster, transition_length, wall_collision_freq = 1e7) Xe_0 = HallThruster.Xenon(0) diff --git a/test/unit_tests/test_misc.jl b/test/unit_tests/test_misc.jl index d75e2760..b1d37b9f 100644 --- a/test/unit_tests/test_misc.jl +++ b/test/unit_tests/test_misc.jl @@ -128,36 +128,3 @@ end @test Aϵ isa SparseArrays.Tridiagonal end - -@testset "Transition functions" begin - - step = HallThruster.StepFunction() - - cutoff = 12 - y1 = 17 - y2 = 101 - - transition_length = 1.0 - offset = 0.0 - - @test step(0, cutoff, y1, y2) == y1 - @test step(cutoff * 2, cutoff, y1, y2) == y2 - - smooth = HallThruster.SmoothIf(transition_length) - - @test smooth(0, cutoff, y1, y2) ≈ y1 - @test smooth(cutoff * 2, cutoff, y1, y2) ≈ y2 - - linear = HallThruster.LinearTransition(transition_length, offset) - - @test linear(0, cutoff, y1, y2) == y1 - @test linear(cutoff * 2, cutoff, y1, y2) == y2 - @test y1 < linear(cutoff + transition_length/10, cutoff, y1, y2) < y2 - @test y1 < linear(cutoff - transition_length/10, cutoff, y1, y2) < y2 - - quadratic = HallThruster.QuadraticTransition(1.0, 0.0) - @test quadratic(0, cutoff, y1, y2) ≈ 17 - @test quadratic(cutoff * 2000, cutoff, y1, y2) ≈ 101 - @test quadratic(cutoff, cutoff, y1, y2) == 17 - @test quadratic(cutoff + transition_length / 10, cutoff, y1, y2) > 1 -end diff --git a/test/unit_tests/test_walls.jl b/test/unit_tests/test_walls.jl index 62fe5688..7f102644 100644 --- a/test/unit_tests/test_walls.jl +++ b/test/unit_tests/test_walls.jl @@ -1,6 +1,6 @@ @testset "Wall loss tests" begin no_losses = HallThruster.NoWallLosses() - + h = HallThruster.edge_to_center_density_ratio() mi = HallThruster.Xenon.m me = HallThruster.me @@ -25,10 +25,10 @@ radial_loss_frequency = [0.0, 0.0, 0.0, 0.0], νew_momentum = [0.0, 0.0, 0.0, 0.0], ) - transition_function = HallThruster.LinearTransition(0.2 * L_ch, 0.0) + transition_length = 0.2 * L_ch config = (; thruster = HallThruster.SPT_100, - transition_function, + transition_length, propellant = HallThruster.Xenon, ncharge = 1, electron_plume_loss_scale = 0.0 ) @@ -98,10 +98,10 @@ @test HallThruster.freq_electron_wall(no_losses, params, 3) == 0.0 @test HallThruster.freq_electron_wall(landmark_losses, params, 2) == 1e7 - @test HallThruster.freq_electron_wall(landmark_losses, params, 3) * params.config.transition_function(params.z_cell[3], params.L_ch, 1.0, 0.0) == 0.0e7 + @test HallThruster.freq_electron_wall(landmark_losses, params, 3) * HallThruster.linear_transition(params.z_cell[3], params.L_ch, params.config.transition_length, 1.0, 0.0) == 0.0e7 γ = HallThruster.SEE_yield(BN, Tev, γmax) - νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr + νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr * h νew = νiw / (1 - γ) params.cache.γ_SEE .= γ @@ -118,17 +118,18 @@ @test Iiw ≈ νiw * HallThruster.e * V_cell * ne @test Iew ≈ νew * HallThruster.e * V_cell * ne - @test HallThruster.freq_electron_wall(sheath_model, params, 2) * params.config.transition_function(params.z_cell[2], L_ch, 1.0, 0.0) ≈ νew - @test HallThruster.freq_electron_wall(sheath_model, params, 3) * params.config.transition_function(params.z_cell[3], L_ch, 1.0, 0.0) ≈ 0.0 + @test HallThruster.freq_electron_wall(sheath_model, params, 2) * HallThruster.linear_transition(params.z_cell[2], L_ch, params.config.transition_length, 1.0, 0.0) ≈ νew + @test HallThruster.freq_electron_wall(sheath_model, params, 3) * HallThruster.linear_transition(params.z_cell[3], L_ch, params.config.transition_length, 1.0, 0.0) ≈ 0.0 - @test HallThruster.wall_power_loss(sheath_model, params, 2) ≈ νew * (2 * (1 - 0.5 * BN.σ₀) * Tev + (1 - γ) * Vs)/γ + @test HallThruster.wall_power_loss(sheath_model, params, 2) ≈ νew * (2 * Tev + (1 - γ) * Vs) @test HallThruster.wall_power_loss(sheath_model, params, 4) ≈ 0.0 end @testset "Ion wall losses" begin - config = (;thruster = HallThruster.SPT_100, propellant = HallThruster.Krypton, ncharge = 2, transition_function = HallThruster.StepFunction()) + config = (;thruster = HallThruster.SPT_100, propellant = HallThruster.Krypton, ncharge = 2, transition_length = 0.0) geom = HallThruster.SPT_100.geometry Δr = geom.outer_radius - geom.inner_radius + h = HallThruster.edge_to_center_density_ratio() Tn = 300.0 Tev = 3.0 @@ -139,7 +140,7 @@ end γ_SEE_max = 1 - 8.3 * sqrt(HallThruster.me/mi) γ = HallThruster.SEE_yield(HallThruster.BoronNitride, Tev, γ_SEE_max) - νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr + νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr * h νew = νiw * γ / (1 - γ) L_ch = HallThruster.geometry_SPT_100.channel_length @@ -257,7 +258,7 @@ end Δz = z_edge[2] - z_edge[1] V_cell = A_ch * Δz u_bohm = sqrt(HallThruster.e * Tev / mi) - νiw = u_bohm / Δr + νiw = u_bohm / Δr * h @test dU[index.ρi[1], i] ≈ -U[index.ρi[1], i] * νiw @test dU[index.ρi[2], i] ≈ -U[index.ρi[2], i] * νiw * sqrt(2) @@ -284,7 +285,7 @@ end params_wall_sheath = (;base_params..., config = config_wall_sheath) γ = HallThruster.SEE_yield(material, Tev, γ_SEE_max) - νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr + νiw = α * sqrt(HallThruster.e * Tev / mi) / Δr * h νew = νiw * γ / (1 - γ) params_wall_sheath.cache.νew_momentum[1:2] .= νew From 1865d5c2c8254748ee4f6c77319146af257899b9 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Wed, 5 Jun 2024 13:53:17 -0400 Subject: [PATCH 12/14] add neutral ingestion multiplier --- src/simulation/boundaryconditions.jl | 5 +---- test/unit_tests/test_boundary_conditions.jl | 9 +++++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/simulation/boundaryconditions.jl b/src/simulation/boundaryconditions.jl index f9a68b38..11aa824b 100644 --- a/src/simulation/boundaryconditions.jl +++ b/src/simulation/boundaryconditions.jl @@ -29,7 +29,7 @@ function left_boundary_state!(bc_state, U, params) bc_state[index.ρn] = mdot_a / channel_area[1] / un # Add ingested mass flow rate at anode - bc_state[index.ρn] += params.background_neutral_density * params.background_neutral_velocity / un + 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] @@ -61,9 +61,6 @@ 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] -= boundary_flux / un diff --git a/test/unit_tests/test_boundary_conditions.jl b/test/unit_tests/test_boundary_conditions.jl index d80b6675..08be6f63 100644 --- a/test/unit_tests/test_boundary_conditions.jl +++ b/test/unit_tests/test_boundary_conditions.jl @@ -17,6 +17,7 @@ ion_temperature = 300.0, background_pressure = 6e-4, background_neutral_temperature = 150.0, + neutral_ingestion_multiplier = 1.5 ) fluids, fluid_ranges, species, species_range_dict, is_velocity_index = HallThruster.configure_fluids(config) @@ -96,13 +97,17 @@ # when ion velocity at left boundary is less than bohm speed, it should be accelerated # to reach bohm speed HallThruster.left_boundary_state!(U_b, U1, params) - @test U_b[index.ρn] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un + @test U_b[index.ρn] ≈ + HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + + background_neutral_density * background_neutral_velocity / un * config.neutral_ingestion_multiplier @test U_b[index.ρiui[1]] / U_b[index.ρi[1]] == -u_bohm_1 @test U_b[index.ρiui[2]] / U_b[index.ρi[2]] == -u_bohm_2 # when ion velocity at left boundary is greater than bohm speed, ions have Neumann BC HallThruster.left_boundary_state!(U_b, U2, params) - @test U_b[index.ρn] ≈ HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + background_neutral_density * background_neutral_velocity / un + @test U_b[index.ρn] ≈ + HallThruster.inlet_neutral_density(config) - (U_b[index.ρiui[1]] + U_b[index.ρiui[2]]) / config.neutral_velocity + + background_neutral_density * background_neutral_velocity / un * config.neutral_ingestion_multiplier @test U_b[index.ρiui[1]] == U_2[index.ρiui[1]] @test U_b[index.ρiui[2]] == U_2[index.ρiui[2]] @test U_b[index.ρiui[1]] / U_b[index.ρi[1]] == -2 * u_bohm_1 From c78ea8e86ac83bac983a4ec7aa8839ddfbc1d24b Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Wed, 5 Jun 2024 13:56:28 -0400 Subject: [PATCH 13/14] split electron and heavy species initialization --- src/simulation/initialization.jl | 46 +++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 13 deletions(-) diff --git a/src/simulation/initialization.jl b/src/simulation/initialization.jl index 3b560d77..51fd052a 100644 --- a/src/simulation/initialization.jl +++ b/src/simulation/initialization.jl @@ -5,13 +5,13 @@ struct DefaultInitialization <: InitialCondition end initialize!(U, params) = initialize!(U, params, params.config.initial_condition) -function initialize!(U, params, model::InitialCondition) +function initialize!(_, _, model::InitialCondition) throw(ArgumentError("Function HallThruster.initialize!(U, params, model::$(typeof(model)) not yet implemented. For InitialCondition types other than DefaultInitialization(), this must be defined by the user!")) end -function initialize!(U, params, ::DefaultInitialization) +function initialize_heavy_species_default!(U, params) (;z_cell, config, index, cache) = params - (;ncharge, anode_Te, cathode_Te, domain, thruster, propellant, discharge_voltage, anode_mass_flow_rate) = config + (;ncharge, anode_Te, domain, thruster, propellant, discharge_voltage, anode_mass_flow_rate) = config mi = propellant.m L_ch = thruster.geometry.channel_length z0 = domain[1] @@ -44,17 +44,11 @@ function initialize!(U, params, ::DefaultInitialization) # Beam neutral density at outlet ρn_1 = 0.01 * ρn_0 - neutral_function = z -> smooth_if(z-z0, L_ch / 2, ρn_0, ρn_1, L_ch / 6) + neutral_function = z -> smooth_if(z-z0, L_ch / 2, ρn_0, ρn_1, L_ch / 6) + # Electron density number_density_function = z -> sum(Z * ion_density_function(z, Z) / mi for Z in 1:ncharge) - Te_baseline = z -> lerp(z, domain[1], domain[2], anode_Te, cathode_Te) - Te_min = min(anode_Te, cathode_Te) - Te_max = (config.discharge_voltage / 10) - Te_width = L_ch/3 - - energy_function = z -> 3/2 * (Te_baseline(z) + (Te_max - Te_min) * exp(-(((z-z0) - L_ch) / Te_width)^2)) - # Fill the state vector for (i, z) in enumerate(z_cell) @@ -65,8 +59,34 @@ function initialize!(U, params, ::DefaultInitialization) U[index.ρiui[Z], i] = ion_density_function(z, Z) * ion_velocity_function(z, Z) end - cache.nϵ[i] = number_density_function(z) * energy_function(z) + cache.ne[i] = number_density_function(z) end - return U + return nothing +end + +function initialize_electrons_default!(params) + (;z_cell, config, cache) = params + (;anode_Te, cathode_Te, domain, thruster) = config + L_ch = thruster.geometry.channel_length + z0 = domain[1] + # Electron temperature + Te_baseline = z -> lerp(z, domain[1], domain[2], anode_Te, cathode_Te) + Te_min = min(anode_Te, cathode_Te) + Te_max = (config.discharge_voltage / 10) + Te_width = L_ch/3 + + # Gaussian Te profile + energy_function = z -> 3/2 * (Te_baseline(z) + (Te_max - Te_min) * exp(-(((z-z0) - L_ch) / Te_width)^2)) + + for (i, z) in enumerate(z_cell) + cache.nϵ[i] = cache.ne[i] * energy_function(z) + end + + return nothing +end + +function initialize!(U, params, ::DefaultInitialization) + initialize_heavy_species_default!(U, params) + initialize_electrons_default!(params) end From 76ca08f95b89c5b5d21f5a787d83b2872ae724b6 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Wed, 5 Jun 2024 14:23:47 -0400 Subject: [PATCH 14/14] update json input --- src/HallThruster.jl | 1 + src/simulation/json.jl | 144 +++++++++++++++++++++ src/simulation/simulation.jl | 199 ----------------------------- test/unit_tests/input_shifted.json | 84 ++++++------ test/unit_tests/input_twozone.json | 76 +++++------ 5 files changed, 219 insertions(+), 285 deletions(-) create mode 100644 src/simulation/json.jl diff --git a/src/HallThruster.jl b/src/HallThruster.jl index 675e5995..f9d4ca3a 100644 --- a/src/HallThruster.jl +++ b/src/HallThruster.jl @@ -62,6 +62,7 @@ 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") diff --git a/src/simulation/json.jl b/src/simulation/json.jl new file mode 100644 index 00000000..84745810 --- /dev/null +++ b/src/simulation/json.jl @@ -0,0 +1,144 @@ +function get_key(json_content, key, default) + if haskey(json_content, key) + return json_content[key] + else + return default + end +end + +function run_simulation(json_content::JSON3.Object, verbose = true) + pressure_z0 = NaN + pressure_dz = NaN + pressure_pstar = NaN + pressure_alpha = NaN + apply_thrust_divergence_correction = true + adaptive = true + + (; + # Design + channel_length, inner_radius, outer_radius, + magnetic_field_file, magnetically_shielded, + propellant, wall_material, + anode_potential, cathode_potential, + anode_mass_flow_rate, + # Simulation + anom_model, cathode_location_m, max_charge, + num_cells, dt_s, duration_s, num_save, + flux_function, limiter, reconstruct, + ion_wall_losses, electron_ion_collisions, + # Parameters + sheath_loss_coefficient, + ion_temp_K, neutral_temp_K, neutral_velocity_m_s, + cathode_electron_temp_eV, inner_outer_transition_length_m, + background_pressure_Torr, background_temperature_K, + ) = json_content + + # Handle optional keys + + # Thruster name + thruster_name = get_key(json_content, :thruster_name, "Unnamed Thruster") + + # Adaptive timestepping + adaptive = get_key(json_content, :adaptive, true) + + # Anom coefficients + anom_model_coeffs = get_key(json_content, :anom_model_coeffs, [0.00625, 0.0625]) + + # Whether inbuilt divergence model is used to correct thrust + apply_thrust_divergence_correction = get_key(json_content, :apply_thrust_divergence_correction, true) + + # neutral ingestion multiplier + neutral_ingestion_multiplier = get_key(json_content, :neutral_ingestion_multiplier, 1.0) + + # Optional parameters for pressure-dependent models + if anom_model == "ShiftedTwoZone" || anom_model == "ShiftedTwoZoneBohm" || + anom_model == "ShiftedMultiBohm" || anom_model == "ShiftedGaussianBohm" + (;pressure_z0, pressure_dz, pressure_pstar, pressure_alpha) = json_content + end + + geometry = Geometry1D(;channel_length, outer_radius, inner_radius) + + if thruster_name == "SPT-100" + bfield_func = HallThruster.B_field_SPT_100 $ (0.016, channel_length) + else + bfield_data = readdlm(magnetic_field_file, ',') + bfield_func = HallThruster.LinearInterpolation(bfield_data[:, 1], bfield_data[:, 2]) + end + + thruster = HallThruster.Thruster(; + name = thruster_name, + geometry = geometry, + magnetic_field = bfield_func, + shielded = magnetically_shielded + ) + + # assign anomalous transport model + anom_model = if anom_model == "NoAnom" + NoAnom() + elseif anom_model == "Bohm" + Bohm(anom_model_coeffs[1]) + elseif anom_model == "TwoZoneBohm" + TwoZoneBohm((anom_model_coeffs[1], anom_model_coeffs[2])) + elseif anom_model == "MultiLogBohm" + MultiLogBohm(anom_model_coeffs) + elseif anom_model == "ShiftedTwoZone" || anom_model == "ShiftedTwoZoneBohm" + coeff_tuple = (anom_model_coeffs[1], anom_model_coeffs[2]) + ShiftedTwoZoneBohm(coeff_tuple, pressure_z0, pressure_dz, pressure_pstar, pressure_alpha) + elseif anom_model == "ShiftedMultiBohm" + N = length(anom_model_coeffs) + ShiftedMultiBohm( + anom_model_coeffs[1:N÷2], anom_model_coeffs[N÷2 + 1], pressure_z0, pressure_dz, pressure_pstar, pressure_alpha + ) + elseif anom_model == "ShiftedGaussianBohm" + ShiftedGaussianBohm( + anom_model_coeffs[1], anom_model_coeffs[2], anom_model_coeffs[3], anom_model_coeffs[4], + pressure_z0, pressure_dz, pressure_pstar, pressure_alpha + ) + end + + config = HallThruster.Config(; + thruster, + propellant = eval(Symbol(propellant)), + anom_model, + domain = (0.0, cathode_location_m), + discharge_voltage = anode_potential - cathode_potential, + anode_mass_flow_rate = anode_mass_flow_rate, + cathode_potential = cathode_potential, + ncharge = max_charge, + wall_loss_model = WallSheath(eval(Symbol(wall_material)), sheath_loss_coefficient), + ion_wall_losses = ion_wall_losses, + cathode_Te = cathode_electron_temp_eV, + LANDMARK = false, + ion_temperature = ion_temp_K, + neutral_temperature = neutral_temp_K, + neutral_velocity = neutral_velocity_m_s, + electron_ion_collisions = electron_ion_collisions, + min_electron_temperature = cathode_electron_temp_eV, + transition_length = inner_outer_transition_length_m, + scheme = HyperbolicScheme(; + flux_function = eval(Symbol(flux_function)), limiter = eval(Symbol(limiter)), reconstruct + ), + background_pressure = background_pressure_Torr * u"Torr", + background_neutral_temperature = background_temperature_K * u"K", + neutral_ingestion_multiplier, + apply_thrust_divergence_correction, + ) + + solution = run_simulation( + config; ncells = num_cells, nsave = num_save, + duration = duration_s, dt = dt_s, verbose = verbose, adaptive + ) + + return solution +end + +function run_simulation(json_path::String; is_path = true, kwargs...) + + if is_path + json_content = JSON3.read(read(json_path, String)) + else + json_content = JSON3.read(json_path) + end + + run_simulation(json_content; kwargs...) +end diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index 92262497..42d1d5a5 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -263,202 +263,3 @@ function run_simulation( return sol end - -function run_simulation(json_content::JSON3.Object; single_section = false, nonstandard_keys = false, verbose = true, adaptive = false) - - pressure_z0 = NaN - pressure_dz = NaN - pressure_pstar = NaN - pressure_alpha = NaN - apply_thrust_divergence_correction = true - - if single_section - if nonstandard_keys - (; - # Design - thruster_name, - channel_length, inner_radius, outer_radius, - magnetic_field_file, magnetically_shielded, - propellant_material, wall_material, - anode_potential, cathode_potential, - anode_mass_flow_rate, - # Simulation - anom_model, cathode_location_m, max_charge, - num_cells, dt_s, duration_s, num_save, - flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, - # Parameters - sheath_loss_coefficient, - ion_temp_K, neutral_temp_K, neutral_velocity_m_s, - cathode_electron_temp_eV, inner_outer_transition_length_m, - background_pressure_Torr, background_temperature_K, - ) = json_content - - propellant = propellant_material - else - (; - # Design - thruster_name, - channel_length, inner_radius, outer_radius, - magnetic_field_file, magnetically_shielded, - propellant, wall_material, - anode_potential, cathode_potential, - anode_mass_flow_rate, - # Simulation - anom_model, cathode_location_m, max_charge, - num_cells, dt_s, duration_s, num_save, - flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, - # Parameters - sheath_loss_coefficient, - ion_temp_K, neutral_temp_K, neutral_velocity_m_s, - cathode_electron_temp_eV, inner_outer_transition_length_m, - background_pressure_Torr, background_temperature_K, - ) = json_content - end - - # Handle optional keys - if haskey(json_content, :adaptive) - adaptive = json_content.adaptive - end - - # Anom coefficients - if (haskey(json_content, :anom_model_coeffs)) - anom_model_coeffs = json_content.anom_model_coeffs - else - anom_model_coeffs = [json_content.anom_coeff_1, json_content.anom_coeff_2] - end - - if (haskey(json_content, :apply_thrust_divergence_correction)) - apply_thrust_divergence_correction = simulation.apply_thrust_divergence_correction - end - - # Optional parameters for pressure-dependent models - if anom_model == "ShiftedTwoZone" || anom_model == "ShiftedTwoZoneBohm" || - anom_model == "ShiftedMultiBohm" || anom_model == "ShiftedGaussianBohm" - (;pressure_z0, pressure_dz, pressure_pstar, pressure_alpha) = json_content - end - else - (;design, simulation, parameters) = json_content - (; - thruster_name, - channel_length, inner_radius, outer_radius, - magnetic_field_file, magnetically_shielded, - propellant, wall_material, - anode_potential, cathode_potential, - anode_mass_flow_rate, - ) = design - - (; - anom_model, cathode_location_m, max_charge, - num_cells, dt_s, duration_s, num_save, - flux_function, limiter, reconstruct, - ion_wall_losses, electron_ion_collisions, - ) = simulation - - (; - anom_model_coeffs, sheath_loss_coefficient, - ion_temp_K, neutral_temp_K, neutral_velocity_m_s, - cathode_electron_temp_eV, inner_outer_transition_length_m, - background_pressure_Torr, background_temperature_K, - ) = parameters - - # Handle optional keys - if (haskey(simulation, :adaptive)) - adaptive = simulation.adaptive - end - - if (haskey(simulation, :apply_thrust_divergence_correction)) - apply_thrust_divergence_correction = simulation.apply_thrust_divergence_correction - end - - # Optional parameters for ShiftedTwoZoneBohm - if anom_model == "ShiftedTwoZone" || anom_model == "ShiftedTwoZoneBohm" - (;pressure_z0, pressure_dz, pressure_pstar, pressure_alpha) = parameters - end - end - - geometry = Geometry1D(;channel_length, outer_radius, inner_radius) - - if thruster_name == "SPT-100" - bfield_func = HallThruster.B_field_SPT_100 $ (0.016, channel_length) - else - bfield_data = readdlm(magnetic_field_file, ',') - bfield_func = HallThruster.LinearInterpolation(bfield_data[:, 1], bfield_data[:, 2]) - end - - thruster = HallThruster.Thruster(; - name = thruster_name, - geometry = geometry, - magnetic_field = bfield_func, - shielded = magnetically_shielded - ) - - anom_model = if anom_model == "NoAnom" - NoAnom() - elseif anom_model == "Bohm" - Bohm(anom_model_coeffs[1]) - elseif anom_model == "TwoZoneBohm" - TwoZoneBohm((anom_model_coeffs[1], anom_model_coeffs[2])) - elseif anom_model == "MultiLogBohm" - MultiLogBohm(anom_model_coeffs) - elseif anom_model == "ShiftedTwoZone" || anom_model == "ShiftedTwoZoneBohm" - coeff_tuple = (anom_model_coeffs[1], anom_model_coeffs[2]) - ShiftedTwoZoneBohm(coeff_tuple, pressure_z0, pressure_dz, pressure_pstar, pressure_alpha) - elseif anom_model == "ShiftedMultiBohm" - N = length(anom_model_coeffs) - ShiftedMultiBohm( - anom_model_coeffs[1:N÷2], anom_model_coeffs[N÷2 + 1], pressure_z0, pressure_dz, pressure_pstar, pressure_alpha - ) - elseif anom_model == "ShiftedGaussianBohm" - ShiftedGaussianBohm( - anom_model_coeffs[1], anom_model_coeffs[2], anom_model_coeffs[3], anom_model_coeffs[4], - pressure_z0, pressure_dz, pressure_pstar, pressure_alpha - ) - end - - config = HallThruster.Config(; - thruster, - propellant = eval(Symbol(propellant)), - anom_model, - domain = (0.0, cathode_location_m), - discharge_voltage = anode_potential - cathode_potential, - anode_mass_flow_rate = anode_mass_flow_rate, - cathode_potential = cathode_potential, - ncharge = max_charge, - wall_loss_model = WallSheath(eval(Symbol(wall_material)), sheath_loss_coefficient), - ion_wall_losses = ion_wall_losses, - cathode_Te = cathode_electron_temp_eV, - LANDMARK = false, - ion_temperature = ion_temp_K, - neutral_temperature = neutral_temp_K, - neutral_velocity = neutral_velocity_m_s, - electron_ion_collisions = electron_ion_collisions, - min_electron_temperature = cathode_electron_temp_eV, - transition_length = inner_outer_transition_length_m, - scheme = HyperbolicScheme(; - flux_function = eval(Symbol(flux_function)), limiter = eval(Symbol(limiter)), reconstruct - ), - background_pressure = background_pressure_Torr * u"Torr", - background_neutral_temperature = background_temperature_K * u"K", - apply_thrust_divergence_correction - ) - - solution = run_simulation( - config; ncells = num_cells, nsave = num_save, - duration = duration_s, dt = dt_s, verbose = verbose, adaptive - ) - - return solution -end - -function run_simulation(json_path::String; is_path = true, kwargs...) - - if is_path - json_content = JSON3.read(read(json_path, String)) - else - json_content = JSON3.read(json_path) - end - - run_simulation(json_content; kwargs...) -end diff --git a/test/unit_tests/input_shifted.json b/test/unit_tests/input_shifted.json index 21da469a..cce7e2e5 100644 --- a/test/unit_tests/input_shifted.json +++ b/test/unit_tests/input_shifted.json @@ -1,47 +1,41 @@ { - "parameters": { - "neutral_temp_K": 300, - "neutral_velocity_m_s": 400.0, - "ion_temp_K": 1000, - "cathode_electron_temp_eV": 1.4254966095382686, - "sheath_loss_coefficient": 0.1643173140884967, - "inner_outer_transition_length_m": 0.0019936840328598983, - "anom_model_coeffs": [ - 0.00625, 0.0625 - ], - "background_pressure_Torr": 7.187258216996673e-05, - "background_temperature_K": 300, - "pressure_dz": 0.005, - "pressure_z0": -0.003, - "pressure_pstar": 3e-05, - "pressure_alpha": 43.0 - }, - "design": { - "thruster_name": "SPT-100", - "inner_radius": 0.035, - "outer_radius": 0.05, - "channel_length": 0.025, - "magnetic_field_file": "unit_tests/bfield_spt100.csv", - "wall_material": "BoronNitride", - "magnetically_shielded": false, - "anode_potential": 300.0, - "cathode_potential": 20.0, - "anode_mass_flow_rate": 3e-6, - "propellant": "Xenon" - }, - "simulation": { - "num_cells": 200, - "dt_s": 1e-09, - "duration_s": 1e-7, - "num_save": 1000, - "cathode_location_m": 0.08, - "max_charge": 3, - "flux_function": "rusanov", - "limiter": "van_leer", - "reconstruct": true, - "ion_wall_losses": true, - "electron_ion_collisions": true, - "anom_model": "ShiftedTwoZone", - "adaptive": true - } + "neutral_temp_K": 300, + "neutral_velocity_m_s": 400.0, + "ion_temp_K": 1000, + "cathode_electron_temp_eV": 1.4254966095382686, + "sheath_loss_coefficient": 0.1643173140884967, + "inner_outer_transition_length_m": 0.0019936840328598983, + "anom_model_coeffs": [ + 0.00625, 0.0625 + ], + "background_pressure_Torr": 7.187258216996673e-05, + "background_temperature_K": 300, + "pressure_dz": 0.005, + "pressure_z0": -0.003, + "pressure_pstar": 3e-05, + "pressure_alpha": 43.0, + "thruster_name": "SPT-100", + "inner_radius": 0.035, + "outer_radius": 0.05, + "channel_length": 0.025, + "magnetic_field_file": "unit_tests/bfield_spt100.csv", + "wall_material": "BoronNitride", + "magnetically_shielded": false, + "anode_potential": 300.0, + "cathode_potential": 20.0, + "anode_mass_flow_rate": 3e-6, + "propellant": "Xenon", + "num_cells": 200, + "dt_s": 1e-09, + "duration_s": 1e-7, + "num_save": 1000, + "cathode_location_m": 0.08, + "max_charge": 3, + "flux_function": "rusanov", + "limiter": "van_leer", + "reconstruct": true, + "ion_wall_losses": true, + "electron_ion_collisions": true, + "anom_model": "ShiftedTwoZone", + "adaptive": true } \ No newline at end of file diff --git a/test/unit_tests/input_twozone.json b/test/unit_tests/input_twozone.json index 9240402f..61b80828 100644 --- a/test/unit_tests/input_twozone.json +++ b/test/unit_tests/input_twozone.json @@ -1,43 +1,37 @@ { - "parameters": { - "neutral_temp_K": 300, - "neutral_velocity_m_s": 400.0, - "ion_temp_K": 1000, - "cathode_electron_temp_eV": 1.4254966095382686, - "sheath_loss_coefficient": 0.1643173140884967, - "inner_outer_transition_length_m": 0.0019936840328598983, - "anom_model_coeffs": [ - 0.00625, 0.0625 - ], - "background_pressure_Torr": 7.187258216996673e-05, - "background_temperature_K": 300 - }, - "design": { - "thruster_name": "SPT-100", - "inner_radius": 0.035, - "outer_radius": 0.05, - "channel_length": 0.025, - "magnetic_field_file": "unit_tests/bfield_spt100.csv", - "wall_material": "BoronNitride", - "magnetically_shielded": false, - "anode_potential": 300.0, - "cathode_potential": 20.0, - "anode_mass_flow_rate": 3e-6, - "propellant": "Xenon" - }, - "simulation": { - "num_cells": 200, - "dt_s": 1e-09, - "duration_s": 1e-7, - "num_save": 1000, - "cathode_location_m": 0.08, - "max_charge": 3, - "flux_function": "rusanov", - "limiter": "van_leer", - "reconstruct": true, - "ion_wall_losses": true, - "electron_ion_collisions": true, - "anom_model": "TwoZoneBohm", - "adaptive": true - } + "neutral_temp_K": 300, + "neutral_velocity_m_s": 400.0, + "ion_temp_K": 1000, + "cathode_electron_temp_eV": 1.4254966095382686, + "sheath_loss_coefficient": 0.1643173140884967, + "inner_outer_transition_length_m": 0.0019936840328598983, + "anom_model_coeffs": [ + 0.00625, 0.0625 + ], + "background_pressure_Torr": 7.187258216996673e-05, + "background_temperature_K": 300, + "thruster_name": "SPT-100", + "inner_radius": 0.035, + "outer_radius": 0.05, + "channel_length": 0.025, + "magnetic_field_file": "unit_tests/bfield_spt100.csv", + "wall_material": "BoronNitride", + "magnetically_shielded": false, + "anode_potential": 300.0, + "cathode_potential": 20.0, + "anode_mass_flow_rate": 3e-6, + "propellant": "Xenon", + "num_cells": 200, + "dt_s": 1e-09, + "duration_s": 1e-7, + "num_save": 1000, + "cathode_location_m": 0.08, + "max_charge": 3, + "flux_function": "rusanov", + "limiter": "van_leer", + "reconstruct": true, + "ion_wall_losses": true, + "electron_ion_collisions": true, + "anom_model": "TwoZoneBohm", + "adaptive": true } \ No newline at end of file