Stochastic Forcing API

StochasticForcing{T, N, A<:AbstractArray}

Stochastic forcing in Fourier space with white-noise temporal correlation. Supports both CPU and GPU architectures.

Mathematical Properties

The forcing F̂(k,t) satisfies:

• ⟨F̂(k,t)⟩ = 0 (zero mean)
• ⟨F̂(k,t) F̂*(k',t')⟩ = Q̂(k) δ(k-k') δ(t-t')/(dt) (white noise)

Implementation

At each timestep, the forcing is computed as: F̂(k) = √(Q̂(k)) · ξ(k) / √(dt)

where ξ(k) is complex white noise with |ξ| = 1 and random phase.

Fields

• forcing_spectrum::A: √Q̂(k) - square root of power spectrum
• energy_injection_rate::T: Target energy injection rate ε
• k_forcing::T: Central forcing wavenumber k_f
• dk_forcing::T: Forcing bandwidth δ_f
• dt::T: Current timestep (for proper scaling)
• domain_size::NTuple{N,T}: Domain size (Lx, Ly, ...)
• field_size::NTuple{N,Int}: Grid size (Nx, Ny, ...)
• wavenumbers::NTuple{N,Vector{T}}: Wavenumber arrays (kx, ky, ...)
• cached_forcing::AbstractArray{Complex{T},N}: Cached forcing (constant within timestep)
• prevsol::Union{Nothing,AbstractArray{Complex{T},N}}: Previous solution (for Stratonovich work)
• rng::AbstractRNG: Random number generator (default: fresh MersenneTwister per instance for thread/parallel safety)
• random_phases::AbstractArray{T,N}: Pre-allocated random phase buffer (on target architecture)
• last_update_time::T: Time of last forcing update
• spectrum_type::Symbol: Type of forcing spectrum
• enforce_hermitian::Bool: Enforce Hermitian symmetry for real-valued fields
• architecture::AbstractArchitecture: CPU() or GPU() architecture

DeterministicForcing{T, N}

Deterministic (non-random) forcing.

Example

# Sinusoidal forcing
forcing = DeterministicForcing(
    (x, y, t, p) -> p[:A] * sin(p[:k] * x) * cos(p[:ω] * t),
    (64, 64);
    parameters = Dict(:A => 1.0, :k => 4.0, :ω => 1.0)
)

generate_forcing!(forcing::StochasticForcing, t::Real, substep::Int=1)

Generate stochastic forcing at time t. Works on both CPU and GPU.

Key Points

1. Forcing is regenerated only on substep 1 - ensures forcing stays constant within a timestep for IMEX and multi-stage methods.

2. Scaling: F̂(k) = √Q̂(k) · exp(i·φ) / √dt, where φ ∼ Uniform[0, 2π) (generated by _fill_random_phases!)

3. Zero mean: The k=0 mode is always set to zero.

4. Hermitian symmetry (optional): F̂(-k) = F̂(k)* so that ifft produces real fields.

5. GPU support: Random phases are generated on the target architecture when possible, then combined with the spectrum using broadcasted operations.

Arguments

• forcing: StochasticForcing configuration
• t: Current simulation time
• substep: Current substep (1 for first substep)

Returns

The cached forcing array (modified in-place).

generate_forcing!(forcing::StochasticForcing, t::Real)

Generate forcing without substep tracking. Equivalent to substep=1.

generate_forcing!(forcing::DeterministicForcing, grid, t::Real)

Evaluate deterministic forcing at time t on the given grid.

apply_forcing!(rhs::AbstractArray, forcing::StochasticForcing, t::Real, substep::Int=1)

Add stochastic forcing to the RHS in spectral space. Works on both CPU and GPU.

Arguments

• rhs: Right-hand side array (modified in-place)
• forcing: StochasticForcing configuration
• t: Current simulation time
• substep: Current substep number

set_dt!(forcing::StochasticForcing, dt::Real)

Update the timestep for Stratonovich scaling.

Call this when dt changes (e.g., adaptive timestepping).

reset_forcing!(forcing::StochasticForcing)

Reset the forcing cache, causing regeneration on next call. Works on both CPU and GPU.

store_prevsol!(forcing::StochasticForcing, sol::AbstractArray)

Store the current solution for Stratonovich work calculation. Works on both CPU and GPU arrays.

Call this at the beginning of each timestep, before advancing.

work_stratonovich(forcing::StochasticForcing, sol::AbstractArray)

Compute work done by forcing using Stratonovich interpretation. Works on both CPU and GPU arrays.

Formula

W = Re⟨(ψⁿ + ψⁿ⁺¹)/2 · ΔF̂*⟩

where ΔF̂ = √Q̂ · ξ · √dt is the forcing increment over dt.

This correctly accounts for the correlation between forcing and response.

Arguments

• forcing: StochasticForcing with prevsol stored
• sol: Current solution ψⁿ⁺¹

Returns

Work done during this timestep (scalar, units of energy).

work_ito(forcing::StochasticForcing, sol::AbstractArray)

Compute work done by forcing using Itô interpretation. Works on both CPU and GPU arrays.

Formula

W_Itô = Re⟨ψⁿ · ΔF̂*⟩ + ε · dt

where ΔF̂ = √Q̂ · ξ · √dt is the forcing increment.

The drift correction ε · dt accounts for the Itô-Stratonovich conversion. In Itô calculus, ψⁿ is independent of Fⁿ⁺¹, so ⟨ψⁿ · ΔF̂⟩ = 0. The drift ensures ⟨WItô⟩ = ⟨WStratonovich⟩ = ε · dt.

mean_energy_injection_rate(forcing::StochasticForcing)

Return the target (mean) energy injection rate ε.

This is the ensemble average of work done per unit time.

instantaneous_power(forcing::StochasticForcing, sol::AbstractArray)

Compute instantaneous power input P = Re⟨ψ · F̂*⟩ / A. Works on both CPU and GPU arrays.

This is the correlation between the solution and forcing at a given instant. For white noise forcing, the expected value depends on which solution is passed:

• ⟨P⟩ = 0 if sol is the solution BEFORE the forcing was applied (independent)
• ⟨P⟩ = ε if sol is the MIDPOINT (ψⁿ + ψⁿ⁺¹)/2
• ⟨P⟩ = 2ε if sol is the solution AFTER forcing (includes full response)

Note: For the Stratonovich-consistent time-averaged power over the timestep, use work_stratonovich(forcing, sol) / forcing.dt instead.

Returns

Instantaneous power (energy per unit time).

forcing_enstrophy_injection_rate(forcing::StochasticForcing)

Compute mean enstrophy injection rate (for 2D turbulence).

η = ∑_k |k|² Q̂(k) / (2 · domain_area)

Note: This assumes the forcing spectrum Q̂(k) corresponds to direct field forcing. For vorticity forcing, enstrophy injection is simply ε. For streamfunction forcing, this formula gives the enstrophy injection rate.

get_forcing_spectrum(forcing::StochasticForcing)

Return the forcing amplitude spectrum √Q̂(k).

get_cached_forcing(forcing::StochasticForcing)

Return the current cached forcing F̂(k).

get_cached_forcing(state::TimestepperState)

Get the cached forcing array. Returns nothing if no forcing is configured.