Postprocessing

Processors

Processors observe the solution inside solve_unsteady after every time step, without having to store the full time history. They are passed as a named tuple, and their outputs are returned with the same field names:

state, outputs = solve_unsteady(;
    setup, start, tlims, params,
    processors = (;
        log = timelogger(; nupdate = 100),
        ehist = realtimeplotter(; setup, plot = energy_history_plot, nupdate = 10),
        vtk = vtk_writer(; setup, dir = "output"),
    ),
)

The built-in processors are timelogger, realtimeplotter/animator (live plots and animations, require a Makie backend), vtk_writer (for ParaView), fieldsaver, observefield, and observespectrum. Custom processors are created with processor.

Note that the state observable passed to a processor contains fields living on the device; for GPU simulations you may have to move them to the host with Array before processing.

Plotting

Loading a Makie backend (e.g. CairoMakie or GLMakie) activates the plotting extension, which provides fieldplot, energy_history_plot, and energy_spectrum_plot. These work both as standalone plots of a state and as plot arguments to realtimeplotter. Derived field quantities such as vorticity and the Q-criterion are computed with vorticity, qcrit, and kinetic_energy; since the velocity components live on staggered positions, they are interpolated to the pressure points (interpolate_u_p) for visualization.

Fields can also be exported to VTK files with save_vtk and viewed in ParaView.

Energy spectra and turbulence statistics

For uniform periodic grids, energyspectrum computes the kinetic energy spectrum $\hat{e}(k)$, which can be compared with the theoretical Kolmogorov scaling $k^{-5/3}$ in 3D (or $k^{-3}$ in 2D) [15]. turbulence_statistics computes quantities such as the Taylor micro-scale and Kolmogorov length scale. The FFT-based machinery is RFFT-based and accounts for the missing conjugate modes; reuse spectral_stuff when computing many spectra on the same grid.

API

Processors

IncompressibleNavierStokes.animator Function

Animate a plot of the solution every update iteration. The animation is saved to path, which should have one of the following extensions:

• ".mkv"

• ".mp4"

• ".webm"

• ".gif"

The plot is determined by a plotter processor. Additional kwargs are passed to plot.

source

IncompressibleNavierStokes.energy_history_plot Function

Create energy history plot.

source

IncompressibleNavierStokes.energy_spectrum_plot Function

Create energy spectrum plot. The energy at a scalar wavenumber level $\kappa \in \mathbb{N}$ is the sum over the wavenumber shell $\kappa \le \parallel k{\parallel }_2<\kappa +1$ (see energyspectrum).

Keyword arguments:

• sloperange = [0.6, 0.9]: Percentage (between 0 and 1) of x-axis where the slope is plotted.

• slopeoffset = 1.3: How far above the energy spectrum the inertial slope is plotted.

source

IncompressibleNavierStokes.fieldplot Function

Plot state field in pressure points. If state is Observable, then the plot is interactive.

Available fieldnames are:

• 1, 2, or 3: velocity component,

• :velocity,

• :velocitynorm,

• :vorticity (default in 2D),

• :qcrit (default in 3D),

• :temperature.

Available plot types for 2D are:

• heatmap (default),

• image,

• contour,

• contourf.

Available plot types for 3D are:

• contour (default).

The alpha value gets passed to contour in 3D.

source

IncompressibleNavierStokes.fieldsaver Method

fieldsaver(; setup, nupdate)

Create processor that stores the solution and time every nupdate time step.

source

IncompressibleNavierStokes.observefield Method

observefield(state; setup, fieldname, psolver)

Observe field fieldname at pressure points.

source

IncompressibleNavierStokes.observespectrum Method

observespectrum(state; setup)

Observe energy spectrum of state (see energyspectrum). Return (; ehat, κ), where ehat is an observable energy spectrum vector.

source

IncompressibleNavierStokes.processor Function

processor(
    initialize
) -> NamedTuple{(:initialize, :finalize), <:Tuple{Any, IncompressibleNavierStokes.var"#316#317"}}
processor(
    initialize,
    finalize
) -> NamedTuple{(:initialize, :finalize), <:Tuple{Any, Any}}

Process results from time stepping. Before time stepping, the initialize function is called on an observable of the time stepper state, returning initialized. The observable is updated every time step.

After timestepping, the finalize function is called on initialized and the final state.

See the following example:

function initialize(state)
    s = Ref(0)
    println("Let's sum up the time steps")
    on(state) do (; n, t)
        println("The summand is $n, the time is $t")
        s[] = s[] + n
    end
    s
end

finalize(s, state) = println("The final sum (at time t=$(state[].t)) is $(s[])")
p = processor(initialize, finalize)

When solved for 6 time steps from t=0 to t=2 the displayed output is

Let's sum up the time steps
The summand is 0, the time is 0.0
The summand is 1, the time is 0.4
The summand is 2, the time is 0.8
The summand is 3, the time is 1.2
The summand is 4, the time is 1.6
The summand is 5, the time is 2.0
The final sum (at time t=2.0) is 15

source

IncompressibleNavierStokes.realtimeplotter Function

Processor for plotting the solution in real time.

Keyword arguments:

• plot: Plot function.

• nupdate: Show solution every nupdate time step.

• displayfig: Display the figure at the start.

• screen: If nothing, use default display. If GLMakie.screen() multiple plots can be displayed in separate windows like in MATLAB (see also GLMakie.closeall()).

• displayupdates: Display the figure at every update (if using CairoMakie).

• sleeptime: The sleeptime is slept at every update, to give Makie time to update the plot. Set this to nothing to skip sleeping.

Additional kwargs are passed to the plot function.

source

IncompressibleNavierStokes.save_vtk Method

save_vtk(state; setup, filename, kwargs...)

Save fields to vtk file.

The kwargs are passed to snapshotsaver.

source

IncompressibleNavierStokes.snapshotsaver Method

snapshotsaver(state; setup, fieldnames, psolver)

In the case of a 2D setup, the velocity field is saved as a 3D vector with a z-component of zero, as this seems to be preferred by ParaView.

source

IncompressibleNavierStokes.timelogger Method

timelogger(
;
    showiter,
    showt,
    showdt,
    showmax,
    showspeed,
    nupdate
) -> @NamedTuple{initialize::IncompressibleNavierStokes.var"#319#322"{Bool, Bool, Bool, Bool, Bool, Int64}, finalize::IncompressibleNavierStokes.var"#321#324"}

Create processor that logs time step information.

source

IncompressibleNavierStokes.vtk_writer Method

vtk_writer(; setup, nupdate, dir, filename, kwargs...)

Create processor that writes the solution every nupdate time steps to a VTK file. The resulting Paraview data collection file is stored in "$dir/$filename.pvd". The kwargs are passed to snapshotsaver.

source

Spectral quantities

IncompressibleNavierStokes.energyspectrum Method

energyspectrum(
    u,
    setup;
    stuff
) -> NamedTuple{(:κ, :ehat), <:Tuple{Any, Any}}

Compute the energy spectrum of the velocity field u. The energy at the scalar wavenumber κ ∈ ℕ is defined as the sum over the wavenumber shell κ ≤ |k| < κ + 1:

$$\hat{e}(\kappa )=\sum _{\kappa \le \parallel k\parallel <\kappa +1}\frac{|\hat{u}(k){|}^2}{2},$$

such that sum(ehat) is the mean kinetic energy density ⟨u_i u_i⟩ / 2 (if all wavenumbers are contained in the shells, see spectral_stuff).

Return (; κ, ehat). The precomputed stuff can be reused between calls.

source

IncompressibleNavierStokes.fftfreq_int Method

fftfreq_int(n, i) -> Any

Integer wavenumber corresponding to index i in an FFT array of size n.

source

IncompressibleNavierStokes.rfft_size Method

rfft_size(Np) -> Any

Size of the RFFT of an array of interior size Np.

source

IncompressibleNavierStokes.spectral_energy_field! Method

spectral_energy_field!(
    ehat_field,
    u,
    setup;
    plan,
    uin,
    uhat
)

Accumulate the spectral kinetic energy sum_α |û_α|² / 2 into the real RFFT-sized array ehat_field, where û_α = rfft(u_α) / prod(Np) (such that the missing-mode-accounted sum over all wavenumbers is the spatial mean kinetic energy density ⟨u_i u_i⟩ / 2).

source

IncompressibleNavierStokes.spectral_stuff Method

spectral_stuff(
    setup;
    kmax
) -> NamedTuple{(:κ, :inds, :energyinds, :Nhat), <:NTuple{4, Any}}

Precompute wavenumber shell indices for computing energy spectra. Shell j contains the wavenumbers κ[j] ≤ |k| < κ[j] + 1, where κ = 0:kmax (all integer wavenumber bins).

Return (; κ, inds, energyinds, Nhat), where

• inds[j]: linear indices of the stored RFFT modes in shell j,

• energyinds[j]: same, but with the indices of modes whose conjugate counterpart is not stored in the RFFT array repeated, such that sum(abs2, view(uhat, energyinds[j])) is the total squared modulus of the full FFT shell,

• Nhat: size of the RFFT array.

source

IncompressibleNavierStokes.spectralsum Method

spectralsum(f) -> Any

Sum over an RFFT-sized array, counting the missing conjugate modes.

source

IncompressibleNavierStokes.squared_wavenumbers Method

squared_wavenumbers(setup) -> Any

Squared integer wavenumber magnitudes |k|² for the RFFT array of interior size Np, flattened to a vector (one entry per linear index).

source

IncompressibleNavierStokes.squared_wavenumbers_dimensional Method

squared_wavenumbers_dimensional(setup) -> Any

Squared dimensional wavenumber magnitudes |2π k / L|² in an RFFT-sized array on the device.

source

IncompressibleNavierStokes.turbulence_statistics Method

turbulence_statistics(
    u,
    setup,
    viscosity
) -> Union{NamedTuple{(:e, :uavg, :diss, :enstrophy, :palinstrophy, :enstrophy_diss, :omega_rms, :l_kra, :l_omega, :t_eddy, :t_ens, :Re, :kmax_etaK), <:NTuple{13, Any}}, NamedTuple{(:e, :uavg, :diss, :l_int, :l_tay, :l_kol, :t_int, :t_tay, :t_kol, :Re_int, :Re_tay, :kmax_eta), <:NTuple{12, Any}}}

Compute turbulence statistics of the velocity field u. The statistics are computed spectrally (RFFT of the interior velocity components), and require a uniform periodic grid. The returned quantities differ between 2D and 3D, since two-dimensional turbulence has different cascades (see the 2D and 3D methods).

source

IncompressibleNavierStokes.turbulence_statistics Method

turbulence_statistics(
    _::IncompressibleNavierStokes.Dimension{2},
    u,
    setup,
    viscosity
) -> NamedTuple{(:e, :uavg, :diss, :enstrophy, :palinstrophy, :enstrophy_diss, :omega_rms, :l_kra, :l_omega, :t_eddy, :t_ens, :Re, :kmax_etaK), <:NTuple{13, Any}}

Two-dimensional turbulence statistics (Kraichnan-Batchelor-Leith dual cascade). The small scales are set by the forward enstrophy cascade, so the relevant dissipation is the enstrophy-dissipation rate and the relevant length is the Kraichnan scale (not the Kolmogorov scale):

• e: mean kinetic energy density ⟨u_i u_i⟩ / 2

• uavg: per-component RMS velocity u' = sqrt(e)

• diss: energy dissipation rate ϵ = 2 ν Ω = ν Σ_k |k|² |û|²

• enstrophy: enstrophy Ω = ⟨ω²⟩ / 2

• palinstrophy: palinstrophy P = ⟨|∇ω|²⟩ / 2

• enstrophy_diss: enstrophy-dissipation rate η_Ω = 2 ν P = ν Σ_k |k|⁴ |û|²

• omega_rms: RMS vorticity ω' = sqrt(2 Ω)

• l_kra: Kraichnan length scale η_K = (ν³ / η_Ω)^(1/6)

• l_omega: Taylor-type length scale u' / ω'

• t_eddy: large-scale strain time 1 / ω'

• t_ens: enstrophy-cascade time scale η_Ω^(-1/3)

• Re: Reynolds number u'² / (ν ω')

• kmax_etaK: resolution indicator k_max η_K (well-resolved DNS: ≳ 1)

source

IncompressibleNavierStokes.turbulence_statistics Method

turbulence_statistics(
    _::IncompressibleNavierStokes.Dimension{3},
    u,
    setup,
    viscosity
) -> NamedTuple{(:e, :uavg, :diss, :l_int, :l_tay, :l_kol, :t_int, :t_tay, :t_kol, :Re_int, :Re_tay, :kmax_eta), <:NTuple{12, Any}}

Three-dimensional homogeneous isotropic turbulence statistics (Kolmogorov scaling) [15]:

• e: mean kinetic energy density ⟨u_i u_i⟩ / 2

• uavg: per-component RMS velocity u' = sqrt(2 e / 3)

• diss: dissipation rate ϵ = ν ⟨∇u:∇u⟩ = ν Σ_k |k|² |û|²

• l_int: integral scale estimate L = u'³ / ϵ

• l_tay: Taylor microscale λ = sqrt(15 ν / ϵ) u'

• l_kol: Kolmogorov length scale η = (ν³ / ϵ)^(1/4)

• t_int: large-eddy turnover time L / u'

• t_tay: Taylor time scale λ / u'

• t_kol: Kolmogorov time scale τ_η = sqrt(ν / ϵ)

• Re_int: integral-scale Reynolds number L u' / ν

• Re_tay: Taylor-scale Reynolds number λ u' / ν

• kmax_eta: resolution indicator k_max η (well-resolved DNS: ≳ 1.5)

source

Utils

IncompressibleNavierStokes.enzyme_wrap Function

Wrap a function to return nothing, because Enzyme can not handle vector return values.

source

IncompressibleNavierStokes.get_lims Function

get_lims(x) -> Tuple{Any, Any}
get_lims(x, n) -> Tuple{Any, Any}

Get approximate lower and upper limits of a field x based on the mean and standard deviation ($\mu \pm n\sigma$). If x is constant, a margin of 1e-4 is enforced. This is required for contour plotting functions that require a certain range.

source

IncompressibleNavierStokes.get_perminds Method

get_perminds(N) -> Tuple{Any, Any}

Get permutation indices for DCT.

source

IncompressibleNavierStokes.getoffset Method

getoffset(I) -> Any

Get offset from CartesianIndices I.

source

IncompressibleNavierStokes.getval Method

getval(_::Val{x}) -> Any

Get value contained in Val.

source

IncompressibleNavierStokes.plotgrid Function

plotgrid(x, y; kwargs...)
plotgrid(x, y, z)

Plot nonuniform Cartesian grid.

source

IncompressibleNavierStokes.splitseed Method

splitseed(seed, n) -> AbstractArray

Split random number generator seed into n new seeds.

source