Analysis

Analysis tools for computing diagnostics, statistics, and derived quantities.

CFL Condition

Compute stable timesteps based on flow velocity.

The controller wraps the solver, not the problem — it needs the solver's communicator for the global reduction that makes the step agree across ranks.

using Tarang

coords = CartesianCoordinates("x", "y")
dist   = Distributor(coords; dtype=Float64)
xbasis = RealFourier(coords["x"]; size=64, bounds=(0.0, 2π))
ybasis = RealFourier(coords["y"]; size=64, bounds=(0.0, 2π))
domain = Domain(dist, (xbasis, ybasis))

mesh = Tarang.create_meshgrid(domain)
x, y = mesh["x"], mesh["y"]

u = VectorField(domain, "u")
u.components[1]["g"] = @.  sin(x) * cos(y)
u.components[2]["g"] = @. -cos(x) * sin(y)

problem = IVP([u]; namespace=Dict("u" => u))
add_equation!(problem, "∂t(u) = 0")
solver = InitialValueSolver(problem, RK222(); dt=1e-3)

# Create the CFL controller
cfl = CFL(solver;
    initial_dt = 1e-3,  # dt used until the first computation
    safety     = 0.5,   # safety factor (0.3-0.5 typical)
    cadence    = 1,     # recompute every N iterations
    max_change = 1.5,   # largest dt increase per commit
    min_change = 0.5,   # largest dt decrease per commit
    threshold  = 0.1,   # only commit changes larger than 10% (see below)
    max_dt     = 0.1,   # hard upper bound on dt
)

# Register the fields that constrain the step
add_velocity!(cfl, u)

# Compute a timestep…
dt = compute_timestep(cfl)   # 0.0491 = safety * Δx / max|u| here

# …or hand the controller to `run!` and let it do that each iteration
# run!(solver; cfl=cfl, stop_iteration=1000)

Every keyword is also a mutable field, so cfl.max_dt = 0.05 between steps works. There is no min_dt: dt has no floor, and the safety, min_change and threshold knobs are what keep it from collapsing.

threshold is sticky-dt hysteresis. Changing dt invalidates the cached LHS factorization in the implicit solve, so a proposed step within threshold (relative) of the current one is discarded and the current dt reused. Default 0.1; set 0.0 to commit every change.

Advection only by default

With just add_velocity!, the returned dt accounts for advection only. Diffusion the timestepper integrates explicitly — notably an LES eddy viscosity νₑ, whose spatial variation rules out the implicit path — carries its own limit dt ≤ 1/(2 ν_max Σᵢ Δxᵢ⁻²) that nothing enforces here. Register it:

add_diffusivity!(cfl, nu)                             # constant ν
add_diffusivity!(cfl, get_eddy_viscosity(les_model))  # LES νₑ (per-rank slab)

Both limits then apply and the smaller step wins. See Diffusive Limit.

Global Reductions

Compute global statistics across MPI processes.

using MPI

# Create reducer
reducer = GlobalArrayReducer(MPI.COMM_WORLD)

# Maximum value
global_max = reduce_scalar(reducer, local_max, MPI.MAX)

# Sum
global_sum = reduce_scalar(reducer, local_sum, MPI.SUM)

# Mean (requires division by total elements)
global_mean = global_sum / total_elements

Flow Statistics

Kinetic Energy

function compute_kinetic_energy(u, reducer)
    local_energy = 0.0

    for component in u.components
        Tarang.ensure_layout!(component, :g)
        local_energy += sum(get_grid_data(component).^2) / 2
    end

    return reduce_scalar(reducer, local_energy, MPI.SUM)
end

Enstrophy

function compute_enstrophy(u, reducer)
    # For 2D: ω = ∂v/∂x - ∂u/∂y
    # Enstrophy = ∫ ω² dV

    ux, uy = u.components[1], u.components[2]

    # Compute vorticity (simplified)
    # ... derivative calculation ...

    return reduce_scalar(reducer, local_enstrophy, MPI.SUM)
end

Reynolds Number

function compute_reynolds_number(u, nu, L, reducer)
    Tarang.ensure_layout!(u.components[1], :g)

    # RMS velocity
    local_u2 = sum(get_grid_data(u.components[1]).^2)
    global_u2 = reduce_scalar(reducer, local_u2, MPI.SUM)
    u_rms = sqrt(global_u2 / total_points)

    return u_rms * L / nu
end

Heat Transfer

Nusselt Number

function compute_nusselt(T, w, L, kappa, reducer)
    Tarang.ensure_layout!(T, :g)
    Tarang.ensure_layout!(w, :g)

    # Convective heat flux
    local_flux = sum(get_grid_data(T) .* get_grid_data(w))
    global_flux = reduce_scalar(reducer, local_flux, MPI.SUM)

    # Normalize
    flux_mean = global_flux / total_points

    # Nusselt = 1 + convective/conductive
    Nu = 1.0 + flux_mean * L / kappa

    return Nu
end

Spectral Analysis

Energy Spectrum

function compute_spectrum(field, kmax)
    Tarang.ensure_layout!(field, :c)

    # Initialize spectrum bins
    E_k = zeros(kmax)

    # Get wavenumbers
    k = get_wavenumbers(field.bases[1])

    # Bin energy by wavenumber
    for (i, ki) in enumerate(k)
        k_bin = round(Int, abs(ki))
        if 1 <= k_bin <= kmax
            E_k[k_bin] += abs2(get_coeff_data(field)[i])
        end
    end

    return E_k
end

Shell-Averaged 3D Spectrum

function compute_3d_spectrum(u, kmax)
    E_k = zeros(kmax)

    for component in u.components
        Tarang.ensure_layout!(component, :c)

        kx = get_wavenumbers(component.bases[1])
        ky = get_wavenumbers(component.bases[2])
        kz = get_wavenumbers(component.bases[3])

        for i in eachindex(kx), j in eachindex(ky), k in eachindex(kz)
            k_mag = sqrt(kx[i]^2 + ky[j]^2 + kz[k]^2)
            k_bin = round(Int, k_mag)

            if 1 <= k_bin <= kmax
                E_k[k_bin] += abs2(get_coeff_data(component)[i,j,k])
            end
        end
    end

    return E_k
end

Time Series

Recording Diagnostics

# Storage
times = Float64[]
energies = Float64[]
nusselts = Float64[]

# During simulation
while solver.sim_time < t_end
    step!(solver, dt)

    # Record
    push!(times, solver.sim_time)
    push!(energies, compute_kinetic_energy(u, reducer))
    push!(nusselts, compute_nusselt(T, w, L, kappa, reducer))
end

Saving Time Series

if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    using JLD2
    @save "diagnostics.jld2" times energies nusselts
end

Spatial Averages

Horizontal Average

function horizontal_average(field)
    Tarang.ensure_layout!(field, :g)

    # Average over x (first axis)
    mean(get_grid_data(field), dims=1)
end

Volume Average

function volume_average(field, reducer)
    Tarang.ensure_layout!(field, :g)

    local_sum = sum(get_grid_data(field))
    global_sum = reduce_scalar(reducer, local_sum, MPI.SUM)

    return global_sum / total_points
end