Adding an anomalous transport model

Users of HallThruster may define their own models by defining a custom subtype of AnomalousTransportModel and a few methods.

Static transport model

Suppose we want to implement $nu_{AN} = \beta \omega_{ce}$ (classic Bohm diffusion). This is a fixed anomalous transport model and does not change as the simulation progresses. We would first define our type:

using HallThruster

struct BohmDiffusion <: AnomalousTransportModel
    β::Float64
end

We then need to define the function which computes the anomalous transport in each cell. This is a mutating function which takes two arguments: the vector of anomalous collision frequency values to be updated, and the solver params.

function (model::BohmDiffusion)(νan, params, config)
    e = HallThruster.e
    me = HallThruster.me
    B = params.cache.B

    for i in eachindex(νan)
        ωce = e * B[i] / me
        νan[i] = model.β * ωce
    end

    return νan
end

We can now set anom_model = BohmDiffusion in our Config struct.

Time-varying transport model

Bohm diffusion is a very simple transport model that does not change over time. More complex, physics-based models may depend self-consistently on local plasma properties and evaolve in time

Let's compute an anomalous transport that depends on energy density of electrostatic waves in the plasma. This model derives from Lafleur, Chabert, and Balruud [1] and has the following form:

\[\begin{aligned} \nu_{AN} = K \frac{\nabla \cdot \mathbf{u}_i W}{m_e n_e c_s v_{de}} \end{aligned}\]

In this expression, $K$ is a tunable coefficient, $u_i$ is the ion velocity, $W = n_e k_B T_e$ is the wave energy density, $m_e$ is the electron mass, $n_e$ is the electron number density, $c_s$ is the ion sound speed, and $v_{de}$ is the electron azimuthal drift speed. Let's say we want to save the wave energy density in addition to the anomalous collision frequency, for later analysis. We begin by defining the model:

struct LafleurModel <: AnomalousTransportModel
    K::Float64
end

Next, we add a method to the num_anom_variables function. This function is part of the Anomalous transport interface and tells HallThruster how many extra arrays we need to allocate during simulation initializion. These arrays store one value per cell. Since we want to save the wave energy density, we need 1 anomalous transport variable.

num_anom_variables(::LafleurModel) = 1

Now, we define the behavior of the model in a function, just as we did for Bohm. Since the model is based on assuming the wave energy convects with the ions, we will use upwind differencing for the gradient.

function (model::LafleurModel)(νan, params, config)
    (;grid, cache) = params
    z = grid.cell_centers
    mi = config.propellant.m
    K = model.K
    e = HallThruster.e
    me = HallThruster.me
    (;ne, Tev, ue, ui, νe, anom_variables) = cache

    ncells = length(νan)

    for i in 2:ncells-1

        W = e * ne[i] * Tev[i]
        Hall_param = e * B[i] / me / νe[i]
        vde = Hall_param * ue[i]
        cs = sqrt(e * Tev[i] / mi)

        # Upwind differencing of gradient term
        if ui > 0
            dz = z[i] - z[i-1]
            W_left = e * ne[i-1] * Tev[i-1]
            grad_ui_W = (ui[1, i] * W[i] - ui[1, i-1] * W_left) / dz
        else
            dz = z[i+1] - z[i]
            W_right = e * ne[i+1] * Tev[i+1]
            grad_ui_W = (ui[1, i+1] * W_right - ui[i] * W[i]) / dz
        end

        # Save W to cache.anom_variables[1]
        anom_variables[1][i] = W

        # Return anomalous collision frequency
        νan[i] = abs(K * grad_ui_W / (me * cs * vde * ne))
    end

    # Neumann boundary condition for all variables
    anom_variables[1][1] = anom_variables[1][2]
    anom_variables[1][end] = anom_variables[1][end-1]
    νan[1] = νan[2]
    νan[end] = νan[end-1]

    return νan
end

Once the solution has been generated by calling run_simulation with appropriate arguments, we can retrieve the wave energy density at frame i as

W = solution.frames[i].anom_variables[1]

Models like this can become unstable quickly, so HallThruster also supports the ability to smooth the anomalous transport. This can be accomplished using the anom_smoothing_iters key of the Config, which determines the number of times a smoothing filter will be applied to the results of the transport calculation.

Advanced usage: implementing additional PDEs for anomalous transport

The most complex anomalous transport might require the solution of one or more additional PDEs. This can be accomplished using the same machinery we used above. We demonstrate this by solving the scalar advection equation using first-order upwind differencing in space and forward Euler integration in time, with periodic boundary conditions. The scalar advection equation is given by:

\[\begin{aligned} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 \end{aligned}\]

To begin, we define the model struct and the number of variables we need.

using HallThruster, Plots

struct ScalarAdvection{F} <: HallThruster.AnomalousTransportModel
    advection_velocity::Float64  # Advection advection_velocity
    initializer::F  # Initialization function
end

# Save two auxilliary variables
#   1) Advected quantity u
#   2) gradient of advected quantity (du/dz)
HallThruster.num_anom_variables(::ScalarAdvection) = 2

Next, we define the model function:

function (model::ScalarAdvection)(νan, params, config)
    ncells = length(νan)

    # Extract variables from params
    cache = params.cache
    z = params.z_cell
    u = cache.anom_variables[1]
    du_dz = cache.anom_variables[2]
    a = model.advection_velocity
    dt = params.dt

    if params.iteration[] < 1
        # Initialize
        model.initializer(u, z)
    else
        for i in eachindex(νan)
            # Setup for periodic boundary conditions
            if i == 1
                i_minus_1 = ncells
                dz_minus = z[2] - z[1]
            else
                i_minus_1 = i - 1
                dz_minus = z[i] - z[i-1]
            end

            if i == ncells
                i_plus_1 = 1
                dz_plus = z[2] - z[1]
            else
                i_plus_1 = i + 1
                dz_plus = z[i+1] - z[i]
            end

            # Update gradient using upwind differencing
            if a > 0
                du_dz[i] = (u[i] - u[i_minus_1]) / dz_minus
            else
                du_dz[i] = (u[i_plus_1] - u[i]) / dz_plus
            end
        end

        # Update advected quantity
        for i in eachindex(νan)
            u[i] -= a * du_dz[i] * dt
        end
    end

    # Return a two-zone bohm anomalous transport result,
    # since we don't really care about the anomalous transport in this case.
    # In a more general case, you would set `νan` based on the results of the PDE
    return HallThruster.TwoZoneBohm(1/160, 1/16)(νan, params)
end

And that's it! Now all there is to do is define our simulation parameters and run!

# Set up config
advection_velocity = 1e5
L = 0.08

config = HallThruster.Config(
    domain = (0.0, L),
    anom_model = ScalarAdvection(advection_velocity, initializer),
    thruster = HallThruster.SPT_100,
    discharge_voltage = 300.0,
    anode_mass_flow_rate = 5e-6
)

# Define dt such that CFL condition is obeyed for scalar advection
# Note -- this CFL number is different from the CFL number used by `HallThruster`
dx = L / ncells
CFL = 0.9
dt = min(1e-8, dx * CFL / advection_velocity)
nsteps = 1000
duration = nsteps * dt
simparams = SimParams(
    adaptive = false,
    nsave = nsteps,
    grid = HallThruster.EvenGrid(200)
    dt = dt, duration = duration,
)
ncells = 200

# Run simulation
solution = HallThruster.run_simulation(config, simparams)

# Extract variables from solution
z = solution[:z]
u = [frame.anom_variables[1] for frame in solution.frames]

We can now use a plotting package of our choice to visualize the results to make sure everything worked well. In this case we will use Plots.

using Plots

# Time needed to transit the domain
t_transit = L / advection_velocity

# Number of periods
num_periods = floor(Int, nsteps * dt / t_transit)

# Plot results
p = plot(; framestyle = :box, xlabel = "x", ylabel = "u", title = "First order upwind for scalar advection")
for i in 0:num_periods
    index = round(Int, i * t_transit / dt) + 1
    plot!(
        p, z, u[index], label = "After $i periods",
        linecolor = cgrad(:turbo, num_periods + 1, categorical = true)[i+1]
    )
end
display(p)

This looks correct! In this case, we haven't coupled our PDE solution to the anomalous transport, but one could easily do this. In the same way, systems of two, three, or more coupled PDEs can be solved and related to the anomalous collision frequency.