Source terms

struct SourceTerm{F, M, T, V}

SourceTerm(f, mode = Unsplit(); timescale = nothing)

A source term for the conserved equations, written as one plain function of position, time, and the whole cell — no kernels, no protocol methods. Pass it (or the bare function) to Simulation(...; sources):

f(x1, x2, x3, t, w) -> (member = (slot = rate, …), …)

w is a nested named snapshot of every member's primitives at the cell center: w.gas holds (ρ, v1, v2, v3), plus e for the ideal equations of state and B1, B2, B3 for MHD; w.dust (when the simulation carries dust) is an N-tuple of (ρ, v1, v2, v3), one entry per species, indexed as w.dust[s].

The return names only what the source touches — every unnamed member and slot is zero. Slots are conserved-variable rates, evaluated from primitives: gas takes ρ, m1, m2, m3 and E (total energy, ideal equations of state only); dust takes an N-tuple of partial (ρ, m1, m2, m3) named tuples aligned with w.dust; dust2 = (ρ = rate,) is single-species sugar for species 2 (mixing dust and dust<s> in one return is rejected). Magnetic slots do not exist: field injection goes through EMFSource, which keeps $∇·B = 0$. Every return is validated by a construction-time probe; the kernels trust.

heating(x1, x2, x3, t, w) = (gas = (E = w.gas.ρ / 2,),)

function gravity(x1, x2, x3, t, w)                    # every fluid member
    gas  = (m3 = -w.gas.ρ * g, E = -w.gas.ρ * w.gas.v3 * g)
    # N-generic; ntuple(Val), not tuple map — map's >32-entry fallback breaks GPU
    # compilation
    dust = ntuple(s -> (m3 = -w.dust[s].ρ * g,), Val(length(w.dust)))
    return (; gas, dust)
end

condensation(x1, x2, x3, t, w) =                      # conservative pair
    (gas = (ρ = -k * w.gas.ρ,), dust1 = (ρ = k * w.gas.ρ,))

Keep one return shape per function — the kernels compile f in, and a source branching to differently-shaped returns is type-unstable on the GPU. The idiom for switching a source on and off is the indicator multiply, S * (x1 < 0). Anything richer than a pointwise rate of existing variables — sub-cycling, per-cell memory, workspace, neighbor access, new variables — is a component on the TimeSteppers protocol, one rung up.

mode picks the injection point — Unsplit (every stage, stepper order), Split (once per cycle, after the stages), or Driving (once per cycle, before the stages), the latter two with an optional super-cycling interval. Split and driving updates restore the interior primitives of the members they touch; unsplit sources ride the stage's refresh. Time-dependent sources see the cycle-start time t (t + Δt in split updates) — stage abscissae are not tracked.

timescale, if given, has the same (x1, x2, x3, t, w) signature, returns the source's local timescale, and limits the timestep to its global minimum — fold any safety fraction into the function itself (e.g. 0.1 * w.gas.e / abs(ė) for a cooling time).

On the GPU, f and timescale compile into the update kernels: anything they capture must be a plain value (or a device array) — not an untyped global.

using Aether

heating(x1, x2, x3, t, w) = (gas = (E = w.gas.ρ / 2,),)
SourceTerm(heating, Split())

# output
SourceTerm(heating, Split())

struct Unsplit

Unsplit()

Integrate a SourceTerm at every explicit stage with weight β[stage] * Δt and the stepper's order. This is the default mode.

struct Split{T}

Split(; interval = nothing)

Apply one forward-Euler u += Δt S update after each cycle's stage loop. With interval = T, apply only when crossing a multiple of T, weighted by T (⌊(t + Δt)/T⌋ - ⌊t/T⌋).

struct Driving{T}

Driving(; interval = nothing)

Injection mode of a SourceTerm: one forward-Euler update u += Δt S over the full timestep, applied once per cycle before the stage loop. interval super-cycles the kick exactly as for Split.

UniformGravity(
    g1,
    g2,
    g3
) -> SourceTerm{F, Unsplit, Nothing, Nothing} where F<:Aether.Source.UniformGravityRate

Constant gravitational acceleration g⃗ = (g1, g2, g3) on every fluid member — gravity is universal, there is no per-member opt-out — as an unsplit SourceTerm for Simulation(...; sources): the gas momenta gain ρ g⃗ plus, for the ideal equations of state, the work rate ρ v⃗·g⃗ in the total energy; each dust species' momenta gain ρₛ g⃗. Build it with the mesh's float type.

using Aether

UniformGravity(0.0, 0.0, -0.1)

# output
SourceTerm(UniformGravityRate(0.0, 0.0, -0.1), Unsplit())

ISMCooling(
    mesh::Mesh{FT},
    eos::Union{IdealHydro, IdealHydroS, IdealMHD, IdealMHDS},
    units::CodeUnits;
    heating_rate,
    safety,
    mode
) -> Union{SourceTerm{F, Unsplit, T, Nothing} where {F<:(Aether.Source.ISMCoolingRate{_A, _B, Val{false}} where {_A, _B}), T<:(Aether.Source.ISMCoolingTime{_A, R} where {_A, R<:(Aether.Source.ISMCoolingRate{_A, _B, Val{false}} where {_A, _B})})}, SourceTerm{F, Unsplit, T, Nothing} where {F<:(Aether.Source.ISMCoolingRate{_A, _B, Val{true}} where {_A, _B}), T<:(Aether.Source.ISMCoolingTime{_A, R} where {_A, R<:(Aether.Source.ISMCoolingRate{_A, _B, Val{true}} where {_A, _B})})}}

Optically thin ISM cooling and heating, $\dot e = -n^2 \Lambda(T) + n\,\Gamma$, as a SourceTerm, with ism_cooling_rate as $\Lambda$ (the Koyama & Inutsuka 2002 fit below $10^{4.2}$ K, SPEX and CGOLS above) and constant photoelectric heating $\Gamma =$ heating_rate erg s⁻¹ (the KI02 value by default). Its timescale limits the timestep to safety local cooling/heating times, $\tau =$ safety $\, e/|\dot e|$.

The gas state is scaled to Kelvin and cm⁻³ through units, a CodeUnits — choose its mean_molecular_weight for the gas phase being modeled (≈ 1.27 for neutral atomic ISM). mode picks the injection point: Unsplit per stage or Split once per cycle.

using Aether

mesh = Mesh(CPU(); size = (8, 8, 8), extent = (1, 1, 1))
eos = IdealHydro(5/3, 1e-12, 1e-10)
units = CodeUnits(length = Units.pc, velocity = 1e5, density = 1.27 * Units.mᵤ,
                  mean_molecular_weight = 1.27)
ISMCooling(mesh, eos, units)

# output
SourceTerm(ISMCoolingRate(KI02 + SPEX + CGOLS), Unsplit(); timescale = 0.3 × cooling time)

struct ImplicitISMCooling{FT, V}

ImplicitISMCooling(mesh, eos, units; heating_rate = 2.0e-26,
                   temperature_ceiling = Inf, tolerance = 1.0e-6,
                   maximum_iterations = 50)

Optically thin ISM cooling and heating, $\dot e = -n^2\Lambda(T) + n\Gamma$, solved implicitly: a per-cell backward-Euler step over the full clock.dt, applied through split_update! once per cycle after the stage loop. Unconditionally stable, so — unlike ISMCooling, whose explicit treatment caps dt at a fraction of the local cooling time and stalls in cold gas — it imposes no timestep constraint: stiff cells relax to their thermal-equilibrium temperature in one step, non-stiff cells integrate at first order, and an explicit stepper (RK2/RK3) marches at the hyperbolic CFL limit with no IMEX and no sub-cycling.

heating_rate is the constant photoelectric $\Gamma$ [erg s⁻¹]; the cooling curve, units scaling, and equilibrium match ISMCooling. The internal-energy floor is eos's pressure floor / (γ − 1). temperature_ceiling is an upper cap in Kelvin (default Inf, no cap): a density-dependent pressure cap $T \le$ temperature_ceiling that discards the unphysical hot gas left by strong turbulent shocks in a warm-ISM model, and bounds the maximum sound speed (so the CFL timestep does not collapse in shock-heated cells). tolerance and maximum_iterations bound the safeguarded Newton iteration.

using Aether

mesh = Mesh(CPU(); size = (8, 8, 8), extent = (1, 1, 1))
eos = IdealHydro(5/3, 1e-12, 1e-10)
units = CodeUnits(length = Units.pc, velocity = 1e5, density = 1.27 * Units.mᵤ,
                  mean_molecular_weight = 1.27)
ImplicitISMCooling(mesh, eos, units; temperature_ceiling = 2e4)

# output
ImplicitISMCooling(backward-Euler, KI02 + SPEX + CGOLS)

struct TurbulenceDriving{FT, H, D}

TurbulenceDriving(mesh; dedt, correlation_time = 0, drive_interval = 0,
                  solenoidal_weight = 1, nlow = 1, nhigh = 2,
                  nmid = (nlow + nhigh) / 2, profile = 0, seed = 1)

Random large-scale forcing driving turbulence at the energy injection rate dedt, applied through the drive! slot once per cycle before the stage loop. The acceleration field is a sum over the discrete Fourier modes of the periodic box,

$F(x) = \sum_n A_n \left[ a_n \cos(k_n \cdot x) - b_n \sin(k_n \cdot x) \right],$

whose phase vectors $a_n, b_n$ evolve by an Ornstein–Uhlenbeck process with correlation time correlation_time (0 = white noise, a fresh pattern every kick) and are Helmholtz-projected to $\zeta P_\perp + (1-\zeta) P_\parallel$ with $\zeta =$ solenoidal_weight — 1 is divergence-free driving, 0 purely compressive. Each kick then

1. removes the mean acceleration, so no net momentum is ever injected, and
2. rescales the field so the box-integrated kinetic-energy input over the kick is exactly dedt × volume × Δt,

before applying $\Delta v = F\,\Delta t$ with its exact kinetic energy increment $\Delta E = m \cdot \Delta v + \rho |\Delta v|^2/2$.

The driven band holds wavenumbers nlow $\le |n| \le$ nhigh (in units of $2\pi/L$), weighted by amplitude profile: 0 — the parabola $1 - 4(k - k_c)^2/k_m^2$ of Schmidt et al. (2006) centered on nmid (vanishing at the band edges); 1 — $1 - (k - k_{low})^2/\delta k^2$; 2 — $1 - (k - \delta k)^2/(2\delta k)^2$. Modes with non-positive weight are dropped.

drive_interval > 0 switches to impulsive driving: kicks fire only when the clock passes multiples of the interval, each integrating over the full interval. Phases update on the host from Xoshiro(seed) — identical on every rank, so the force field is MPI-consistent; the two global reductions go through global_sum. Two- and three-dimensional meshes only.

using Aether

mesh = Mesh(CPU(); size = (16, 16, 16), extent = (1, 1, 1))
TurbulenceDriving(mesh; dedt = 1.0)

# output
TurbulenceDriving(dedt = 1.0, 16 modes)

struct EMFSource{F}

EMFSource(f)

A user-defined electric field injected through the constrained-transport update: f is evaluated at edge centers as f(x1, x2, x3, t) and returns (E1, E2, E3), which is added to the corner electric fields every stage — so the field it drives obeys $\partial_t B \mathrel{-}= \nabla \times E$ and stays divergence-free to round-off by construction. Positions along singleton directions evaluate at the domain center, keeping the two edge layers of a 2D/1D run equal so no spurious cross-plane gradients arise.

Only the constrained-transport (MHD) pipeline calls this slot; in a hydrodynamic simulation the component is inert. On the GPU, f compiles into the injection kernels: anything it captures must be a plain value or a device array.

using Aether

driven(x1, x2, x3, t) = (0, 0, x1)
EMFSource(driven)

# output
EMFSource(driven)