hallmd.models

All models are specified as callable functions in hallmd.models. Currently supported models are based on a three-component feedforward system for a Hall thruster:

1. Cathode - Accounts for interactions of the cathode plasma with the main discharge.
2. Thruster - The primary simulation of the Hall thruster channel discharge and near-field.
3. Plume - Models the far-field expansion of the plasma plume in the vacuum chamber.

Fig 1. The three-component feedforward Hall thruster model (Eckels et al 2024).

Examples of integrated predictive engineering models (PEM) are included in the scripts folder.

cathode_coupling(inputs)

Computes cathode coupling voltage dependence on background pressure.

PARAMETER	DESCRIPTION
inputs	input arrays - P_b, V_a, T_e, V_vac, Pstar, P_T for background pressure (Torr), discharge voltage (V), electron temperature (eV), vacuum coupling voltage (V), and model parameters P* (Torr) and P_T (Torr). TYPE: Dataset

RETURNS	DESCRIPTION
Dataset	output arrays - V_cc for cathode coupling voltage (V).

Source code in src/hallmd/models/cathode.py

def cathode_coupling(inputs: Dataset) -> Dataset:
    """Computes cathode coupling voltage dependence on background pressure.

    :param inputs: input arrays - `P_b`, `V_a`, `T_e`, `V_vac`, `Pstar`, `P_T` for background pressure (Torr), discharge
                   voltage (V), electron temperature (eV), vacuum coupling voltage (V), and model parameters P* (Torr)
                   and P_T (Torr).
    :returns outputs: output arrays - `V_cc` for cathode coupling voltage (V).
    """
    # Load inputs
    PB = inputs['P_b'] * TORR_2_PA          # Background Pressure (Torr)
    Va = inputs['V_a']                      # Anode voltage (V)
    Te = inputs['T_e']                      # Electron temperature at the cathode (eV)
    V_vac = inputs['V_vac']                 # Vacuum coupling voltage model parameter (V)
    Pstar = inputs['Pstar'] * TORR_2_PA     # Model parameter P* (Torr)
    PT = inputs['P_T'] * TORR_2_PA          # Model parameter P_T (Torr)

    # Compute cathode coupling voltage
    V_cc = np.atleast_1d(V_vac + Te * np.log(1 + PB / PT) - (Te / (PT + Pstar)) * PB)
    V_cc[V_cc < 0] = 0
    ind = np.where(V_cc > Va)
    V_cc[ind] = np.atleast_1d(Va)[ind]
    return {'V_cc': V_cc}

current_density(inputs, sweep_radius=1.0)

Compute the semi-empirical ion current density (\(j_{ion}\)) plume model over a 90 deg sweep, with 0 deg at thruster centerline. Also compute the plume divergence angle. Will return the ion current density at 91 points, from 0 to 90 deg in 1 deg increments. The angular locations are returned as j_ion_coords in radians.

PARAMETER	DESCRIPTION
inputs	input arrays - P_b, c0, c1, c2, c3, c4, c5, sigma_cex, I_B0 for background pressure (Torr), plume fit coefficients, charge-exchange cross-section (\(m^2\)), and total initial ion beam current (A). If T is provided, then also compute corrected thrust using the divergence angle. TYPE: Dataset
sweep_radius	the location(s) at which to compute the ion current density 90 deg sweep, in units of radial distance (m) from the thruster exit plane. If multiple locations are provided, then the returned \(j_{ion}\) array's last dimension will match the length of sweep_radius. Defaults to 1 meter. TYPE: float | list DEFAULT: 1.0

RETURNS	DESCRIPTION
	output arrays - j_ion for ion current density (\(A/m^2\)) at the j_ion_coords locations, and div_angle in radians for the divergence angle of the plume. Optionally, T_c for corrected thrust (N) if T is provided in the inputs.

Source code in src/hallmd/models/plume.py

def current_density(inputs: Dataset, sweep_radius: float | list = 1.0):
    """Compute the semi-empirical ion current density ($j_{ion}$) plume model over a 90 deg sweep, with 0 deg at
    thruster centerline. Also compute the plume divergence angle. Will return the ion current density at 91 points,
    from 0 to 90 deg in 1 deg increments. The angular locations are returned as `j_ion_coords` in radians.

    :param inputs: input arrays - `P_b`, `c0`, `c1`, `c2`, `c3`, `c4`, `c5`, `sigma_cex`, `I_B0` for background
                   pressure (Torr), plume fit coefficients, charge-exchange cross-section ($m^2$),
                   and total initial ion beam current (A). If `T` is provided, then
                   also compute corrected thrust using the divergence angle.
    :param sweep_radius: the location(s) at which to compute the ion current density 90 deg sweep, in units of radial
                         distance (m) from the thruster exit plane. If multiple locations are provided, then the
                         returned $j_{ion}$ array's last dimension will match the length of `sweep_radius`. Defaults to
                         1 meter.
    :returns outputs: output arrays - `j_ion` for ion current density ($A/m^2$) at the `j_ion_coords` locations,
                                       and `div_angle` in radians for the divergence angle of the plume. Optionally,
                                       `T_c` for corrected thrust (N) if `T` is provided in the inputs.
    """
    # Load plume inputs
    P_B = inputs['P_b'] * TORR_2_PA  # Background pressure (Torr)
    c0 = inputs['c0']  # Fit coefficients (-)
    c1 = inputs['c1']  # (-)
    c2 = inputs['c2']  # (rad/Pa)
    c3 = inputs['c3']  # (rad)
    c4 = inputs['c4']  # (m^-3/Pa)
    c5 = inputs['c5']  # (m^-3)
    sigma_cex = inputs['sigma_cex']  # Charge-exchange cross-section (m^2)
    I_B0 = inputs['I_B0']  # Total initial ion beam current (A)
    thrust = inputs.get('T', None)  # Thrust (N)
    sweep_radius = np.atleast_1d(sweep_radius)

    # 90 deg angle sweep for ion current density
    alpha_rad = np.linspace(0, np.pi / 2, 91)

    # Neutral density
    n = c4 * P_B + c5  # m^-3

    # Divergence angles
    alpha1 = np.atleast_1d(c2 * P_B + c3)  # Main beam divergence (rad)
    alpha1[alpha1 > np.pi / 2] = np.pi / 2
    alpha2 = alpha1 / c1  # Scattered beam divergence (rad)

    with np.errstate(invalid='ignore', divide='ignore'):
        A1 = (1 - c0) / (
            (np.pi ** (3 / 2))
            / 2
            * alpha1
            * np.exp(-((alpha1 / 2) ** 2))
            * (
                2 * erfi(alpha1 / 2)
                + erfi((np.pi * 1j - (alpha1**2)) / (2 * alpha1))
                - erfi((np.pi * 1j + (alpha1**2)) / (2 * alpha1))
            )
        )
        A2 = c0 / (
            (np.pi ** (3 / 2))
            / 2
            * alpha2
            * np.exp(-((alpha2 / 2) ** 2))
            * (
                2 * erfi(alpha2 / 2)
                + erfi((np.pi * 1j - (alpha2**2)) / (2 * alpha2))
                - erfi((np.pi * 1j + (alpha2**2)) / (2 * alpha2))
            )
        )
        # Broadcast over angles and radii (..., a, r)
        A1 = np.expand_dims(A1, axis=(-1, -2))  # (..., 1, 1)
        A2 = np.expand_dims(A2, axis=(-1, -2))
        alpha1 = np.expand_dims(alpha1, axis=(-1, -2))
        alpha2 = np.expand_dims(alpha2, axis=(-1, -2))
        I_B0 = np.expand_dims(I_B0, axis=(-1, -2))
        n = np.expand_dims(n, axis=(-1, -2))
        sigma_cex = np.expand_dims(sigma_cex, axis=(-1, -2))

        decay = np.exp(-sweep_radius * n * sigma_cex)  # (..., 1, r)
        j_cex = I_B0 * (1 - decay) / (2 * np.pi * sweep_radius**2)

        base_density = I_B0 * decay / sweep_radius**2
        j_beam = base_density * A1 * np.exp(-((alpha_rad[..., np.newaxis] / alpha1) ** 2))
        j_scat = base_density * A2 * np.exp(-((alpha_rad[..., np.newaxis] / alpha2) ** 2))

        j_ion = j_beam + j_scat + j_cex  # (..., 91, r) the current density 1d profile at r radial locations

    # Set j~0 where alpha1 < 0 (invalid cases)
    invalid_idx = np.logical_or(np.any(alpha1 <= 0, axis=(-1, -2)), np.any(j_ion <= 0, axis=(-1, -2)))
    j_ion[invalid_idx, ...] = 1e-20
    j_cex[invalid_idx, ...] = 1e-20

    if np.any(abs(j_ion.imag) > 0):
        LOGGER.warning('Predicted beam current has non-zero imaginary component.')
    j_ion = j_ion.real

    # Calculate divergence angle from https://aip.scitation.org/doi/10.1063/5.0066849
    # Requires alpha = [0, ..., 90] deg, from thruster exit-plane to thruster centerline (need to flip)
    # do j_beam + j_scat instead of j_ion - j_cex to avoid catastrophic loss of precision when
    # j_beam and j_scat << j_cex
    j_non_cex = np.flip((j_beam + j_scat).real, axis=-2)
    den_integrand = j_non_cex * np.cos(alpha_rad[..., np.newaxis])
    num_integrand = den_integrand * np.sin(alpha_rad[..., np.newaxis])

    with np.errstate(divide='ignore', invalid='ignore'):
        num = simpson(num_integrand, x=alpha_rad, axis=-2)
        den = simpson(den_integrand, x=alpha_rad, axis=-2)
        cos_div = np.atleast_1d(num / den)
        cos_div[cos_div == np.inf] = np.nan

    div_angle = np.arccos(cos_div)  # Divergence angle (rad) - (..., r)

    # Squeeze last dim if only a single radius was passed
    if sweep_radius.shape[0] == 1:
        j_ion = np.squeeze(j_ion, axis=-1)
        div_angle = np.squeeze(div_angle, axis=-1)

    ret = {'j_ion': j_ion, 'div_angle': div_angle}

    if thrust is not None:
        thrust_corrected = np.expand_dims(thrust, axis=-1) * cos_div
        if sweep_radius.shape[0] == 1:
            thrust_corrected = np.squeeze(thrust_corrected, axis=-1)
        ret['T_c'] = thrust_corrected

    # Interpolate to requested angles
    # if j_ion_coords is not None:
    #     # Extend to range (-90, 90) deg
    #     alpha_grid = np.concatenate((-np.flip(alpha_rad)[:-1], alpha_rad))               # (2M-1,)
    #     jion_grid = np.concatenate((np.flip(j_ion, axis=-1)[..., :-1], j_ion), axis=-1)  # (..., 2M-1)
    #
    #     f = interp1d(alpha_grid, jion_grid, axis=-1)
    #     j_ion = f(j_ion_coords)  # (..., num_pts)

    # Broadcast coords to same loop shape as j_ion (all use the same coords -- store in object array)
    last_axis = -1 if sweep_radius.shape[0] == 1 else -2
    j_ion_coords = np.empty(j_ion.shape[:last_axis], dtype=object)
    for index in np.ndindex(j_ion.shape[:last_axis]):
        j_ion_coords[index] = alpha_rad

    ret['j_ion_coords'] = j_ion_coords

    return ret

hallthruster_jl(thruster_inputs=None, thruster='SPT-100', config=None, simulation=None, postprocess=None, model_fidelity=(2, 2), output_path=None, version=HALLTHRUSTER_VERSION_DEFAULT, pem_to_julia='default', fidelity_function='default', julia_script=None, run_kwargs='default', shock_threshold=None)

Run a single HallThruster.jl simulation for a given set of inputs. This function will write a temporary input file to disk, call HallThruster.run_simulation() in Julia, and read the output file back into Python. Will return time-averaged performance metrics and ion velocity for use with the PEM.

Note that the specific inputs and outputs described here can be configured using the pem_to_julia dict.

Required configuration

You must specify a thruster, a domain, a mass flow rate, and a discharge voltage to run the simulation. The thruster must be defined in the hallmd.devices directory or as a dictionary with the required fields. The mass flow rate and discharge voltage are specified in thruster_inputs as mdot_a (kg/s) and V_a (V), respectively. The domain is specified as a list [left_bound, right_bound] in the config dictionary. See the HallThruster.jl docs for more details.

PARAMETER	DESCRIPTION
thruster_inputs	named key-value pairs of thruster inputs: P_b, V_a, mdot_a, T_e, u_n, l_t, a_1, a_2, delta_z, z0, p0, and V_cc for background pressure (Torr), anode voltage, anode mass flow rate (kg/s), electron temperature (eV), neutral velocity (m/s), transition length (m), anomalous transport coefficients, and cathode coupling voltage. Will override the corresponding values in config if provided. TYPE: Dataset DEFAULT: None
thruster	the name of the thruster to simulate (must be importable from hallmd.devices, see load_device), or a dictionary that provides geometry and magnetic field information of the thruster to simulate; see the Hallthruster.jl docs. Will override thruster in config if provided. If None, will defer to config. Defaults to the SPT-100. TYPE: Literal['SPT-100'] | str | dict DEFAULT: 'SPT-100'
config	dictionary of configs for HallThruster.jl, see the Hallthruster.jl docs for options and formatting. TYPE: Optional[dict] DEFAULT: None
simulation	dictionary of simulation parameters for HallThruster.jl TYPE: Optional[dict] DEFAULT: None
postprocess	dictionary of post-processing parameters for Hallthruster.jl TYPE: Optional[dict] DEFAULT: None
model_fidelity	tuple of integers that determine the number of cells and the number of charge states to use via ncells = model_fidelity[0] * 50 + 100 and ncharge = model_fidelity[1] + 1. Will override ncells and ncharge in simulation and config if provided. TYPE: tuple DEFAULT: (2, 2)
output_path	base path to save output files, will write to current directory if not specified TYPE: Optional[str | Path] DEFAULT: None
version	version of HallThruster.jl to use; will search for a global hallthruster_{version} environment in the ~/.julia/environments/ directory. Can also specify a specific git ref (i.e. branch, commit hash, etc.) to use from GitHub. If the hallthruster_{version} environment does not exist, an error will be raised -- you should create this environment first before using it. TYPE: str DEFAULT: HALLTHRUSTER_VERSION_DEFAULT
pem_to_julia	a dict mapping of PEM shorthand variable names to a list of keys that maps into the HallThruster.jl input/output data structure. Defaults to the provided PEM_TO_JULIA dict defined in hallmd.models.thruster. For example, {'P_b': ['config', 'background_pressure']} will set config['background_pressure'] = P_b. If specified, will override and extend the default mapping. TYPE: dict | str DEFAULT: 'default'
fidelity_function	a callable that takes a tuple of integers and returns a dictionary of simulation parameters. Defaults to _default_model_fidelity which sets ncells and ncharge based on the input tuple. The returned simulation parameters must be convertable to Julia via the pem_to_julia mapping. The callable should also take in the current json config dict. TYPE: Callable[[tuple[int, ...]], dict] DEFAULT: 'default'
julia_script	path to a custom Julia script to run. The script should accept the input json file path as a command line argument. Defaults to just calling HallThruster.run_simulation(input_file). TYPE: Optional[str | Path] DEFAULT: None
run_kwargs	additional keyword arguments to pass to subprocess.run when calling the Julia script. Defaults to check=True. TYPE: dict DEFAULT: 'default'
shock_threshold	if provided, an error will be raised if the ion velocity reaches a maximum before this threshold axial location (in m) - used to detect and filter unwanted "shock-like" behavior, for example by providing a threshold of half the domain length. If not provided, then no filtering is performed (default). TYPE: float DEFAULT: None

RETURNS	DESCRIPTION
Dataset	dict of Hallthruster.jl outputs: I_B0, I_d, T, eta_c, eta_m, eta_v, and u_ion for ion beam current (A), discharge current (A), thrust (N), current efficiency, mass efficiency, voltage efficiency, and singly-charged ion velocity profile (m/s), all time-averaged.

Source code in src/hallmd/models/thruster.py

def hallthruster_jl(
    thruster_inputs: Dataset = None,
    thruster: Literal["SPT-100"] | str | dict = "SPT-100",
    config: Optional[dict] = None,
    simulation: Optional[dict] = None,
    postprocess: Optional[dict] = None,
    model_fidelity: tuple = (2, 2),
    output_path: Optional[str | Path] = None,
    version: str = HALLTHRUSTER_VERSION_DEFAULT,
    pem_to_julia: dict | str = "default",
    fidelity_function: Callable[[tuple[int, ...]], dict] = "default",
    julia_script: Optional[str | Path] = None,
    run_kwargs: dict = "default",
    shock_threshold: float = None,
) -> Dataset:
    """Run a single `HallThruster.jl` simulation for a given set of inputs. This function will write a temporary
    input file to disk, call `HallThruster.run_simulation()` in Julia, and read the output file back into Python. Will
    return time-averaged performance metrics and ion velocity for use with the PEM.

    Note that the specific inputs and outputs described here can be configured using the `pem_to_julia` dict.

    !!! Warning "Required configuration"
        You must specify a thruster, a domain, a mass flow rate, and a discharge voltage to run the simulation. The
        thruster must be defined in the `hallmd.devices` directory or as a dictionary with the required fields.
        The mass flow rate and discharge voltage are specified in `thruster_inputs` as `mdot_a` (kg/s) and
        `V_a` (V), respectively. The domain is specified as a list `[left_bound, right_bound]` in the
        `config` dictionary. See the
        [HallThruster.jl docs](https://um-pepl.github.io/HallThruster.jl/dev/reference/config/) for more details.

    :param thruster_inputs: named key-value pairs of thruster inputs: `P_b`, `V_a`, `mdot_a`, `T_e`, `u_n`, `l_t`,
                            `a_1`, `a_2`, `delta_z`, `z0`, `p0`, and `V_cc` for background pressure (Torr), anode
                            voltage, anode mass flow rate (kg/s), electron temperature (eV), neutral velocity (m/s),
                            transition length (m), anomalous transport coefficients, and cathode coupling voltage. Will
                            override the corresponding values in `config` if provided.
    :param thruster: the name of the thruster to simulate (must be importable from `hallmd.devices`, see
                     [`load_device`][hallmd.utils.load_device]), or a dictionary that provides geometry and
                     magnetic field information of the thruster to simulate; see the
                     [Hallthruster.jl docs](https://um-pepl.github.io/HallThruster.jl/dev/tutorials/simulation//run/).
                     Will override `thruster` in `config` if provided. If None, will defer to `config`.
                     Defaults to the SPT-100.
    :param config: dictionary of configs for `HallThruster.jl`, see the
                   [Hallthruster.jl docs](https://um-pepl.github.io/HallThruster.jl/dev/reference/config/) for
                   options and formatting.
    :param simulation: dictionary of simulation parameters for `HallThruster.jl`
    :param postprocess: dictionary of post-processing parameters for `Hallthruster.jl`
    :param model_fidelity: tuple of integers that determine the number of cells and the number of charge states to use
                           via `ncells = model_fidelity[0] * 50 + 100` and `ncharge = model_fidelity[1] + 1`.
                           Will override `ncells` and `ncharge` in `simulation` and `config` if provided.
    :param output_path: base path to save output files, will write to current directory if not specified
    :param version: version of HallThruster.jl to use; will
                    search for a global `hallthruster_{version}` environment in the `~/.julia/environments/` directory.
                    Can also specify a specific git ref (i.e. branch, commit hash, etc.) to use from GitHub. If the
                    `hallthruster_{version}` environment does not exist, an error will be raised -- you should create
                    this environment first before using it.
    :param pem_to_julia: a `dict` mapping of PEM shorthand variable names to a list of keys that maps into the
                         `HallThruster.jl` input/output data structure. Defaults to the provided PEM_TO_JULIA dict
                         defined in [`hallmd.models.thruster`][hallmd.models.thruster]. For example,
                         `{'P_b': ['config', 'background_pressure']}` will set `config['background_pressure'] = P_b`.
                         If specified, will override and extend the default mapping.
    :param fidelity_function: a callable that takes a tuple of integers and returns a dictionary of simulation
                              parameters. Defaults to `_default_model_fidelity` which sets `ncells` and `ncharge` based
                              on the input tuple. The returned simulation parameters must be convertable to Julia via
                              the `pem_to_julia` mapping. The callable should also take in the current json config dict.
    :param julia_script: path to a custom Julia script to run. The script should accept the input json file path as
                         a command line argument. Defaults to just calling `HallThruster.run_simulation(input_file)`.
    :param run_kwargs: additional keyword arguments to pass to `subprocess.run` when calling the Julia script.
                       Defaults to `check=True`.
    :param shock_threshold: if provided, an error will be raised if the ion velocity reaches a maximum before this
                            threshold axial location (in m) - used to detect and filter unwanted "shock-like" behavior,
                            for example by providing a threshold of half the domain length. If not provided, then no
                            filtering is performed (default).
    :returns: `dict` of `Hallthruster.jl` outputs: `I_B0`, `I_d`, `T`, `eta_c`, `eta_m`, `eta_v`, and `u_ion` for ion
              beam current (A), discharge current (A), thrust (N), current efficiency, mass efficiency, voltage
              efficiency, and singly-charged ion velocity profile (m/s), all time-averaged.
    """
    if pem_to_julia is None or pem_to_julia == "default":
        pem_to_julia = copy.deepcopy(PEM_TO_JULIA)
    else:
        tmp = copy.deepcopy(PEM_TO_JULIA)
        tmp.update(pem_to_julia)
        pem_to_julia = tmp

    thruster_inputs = thruster_inputs or {}

    # Format PEM inputs for HallThruster.jl
    json_data = _format_hallthruster_jl_input(
        thruster_inputs,
        thruster=thruster,
        config=config,
        simulation=simulation,
        postprocess=postprocess,
        model_fidelity=model_fidelity,
        output_path=output_path,
        pem_to_julia=pem_to_julia,
        fidelity_function=fidelity_function,
    )
    # Get julia environment
    jl_environment = get_jl_env(version) if version is not None else None

    if run_kwargs is None:
        run_kwargs = {}
    elif run_kwargs == "default":
        run_kwargs = {"check": True}

    # Run Julia
    t1 = time.time()
    sim_results = run_hallthruster_jl(json_data, jl_env=jl_environment, jl_script=julia_script, **run_kwargs)
    t2 = time.time()

    # Format QOIs for PEM
    thruster_outputs = _convert_to_pem(sim_results, pem_to_julia)

    # Raise an exception if thrust or beam current are negative (non-physical cases)
    thrust = thruster_outputs.get("T", 0)
    beam_current = thruster_outputs.get("I_B0", 0)
    if thrust < 0 or beam_current < 0:
        raise ValueError(f"Exception due to non-physical case: thrust={thrust} N, beam current={beam_current} A")

    # Raise an exception for ion velocity that exhibits "shock" behavior
    if shock_threshold is not None:
        z_coords = thruster_outputs.get("u_ion_coords")
        ion_velocity = thruster_outputs.get("u_ion")
        if z_coords is not None and ion_velocity is not None:
            if (z_max := z_coords[np.argmax(ion_velocity)]) < shock_threshold:
                raise ValueError(f"Exception due to shock-like behavior: max ion velocity occurs at z={z_max:.3f} m")

    thruster_outputs["model_cost"] = t2 - t1  # seconds

    if output_path is not None:
        output_file = Path(json_data["postprocess"].get("output_file"))
        thruster_outputs["output_path"] = output_file.relative_to(Path(output_path).resolve()).as_posix()

    return thruster_outputs