Large Eddy Simulation (LES) Models

Tarang.jl provides Smagorinsky and Anisotropic Minimum Dissipation (AMD) subgrid-scale closures.

What is Large Eddy Simulation?

LES resolves large eddies and models the effect of unresolved motions. The filter width $\Delta$, typically set from grid spacing, separates these scales.

Filtering the Navier-Stokes Equations

We apply a spatial filter to decompose velocity into resolved ($\bar{u}$) and subgrid ($u'$) parts:

\[u = \bar{u} + u'\]

Filtering the incompressible Navier-Stokes equations gives:

\[\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} = -\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x_i} + \nu \nabla^2 \bar{u}_i - \frac{\partial \tau_{ij}}{\partial x_j}\]

where the subgrid-scale (SGS) stress tensor appears:

\[\tau_{ij} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j\]

This tensor represents the effect of unresolved turbulent motions on the resolved flow. Since we cannot compute $\tau_{ij}$ directly (it involves unresolved velocities), we must model it.

The Closure Problem

The Eddy Viscosity Hypothesis

Most SGS models use the Boussinesq hypothesis, which assumes the SGS stress is proportional to the resolved strain rate:

\[\tau_{ij} - \frac{1}{3}\tau_{kk}\delta_{ij} = -2 \nu_e \bar{S}_{ij}\]

where:

• \[\nu_e\]

  is the eddy viscosity (to be modeled)

• \[\bar{S}_{ij}\]

  is the resolved strain rate tensor:

\[\bar{S}_{ij} = \frac{1}{2}\left(\frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right)\]

This transforms the filtered momentum equation into:

\[\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} = -\frac{1}{\rho}\frac{\partial \bar{p}^*}{\partial x_i} + \frac{\partial}{\partial x_j}\left[2(\nu + \nu_e)\bar{S}_{ij}\right]\]

where we've absorbed the isotropic part into a modified pressure $\bar{p}^*$.

For constant molecular viscosity and incompressible velocity, the viscous term expands to $(\nu + \nu_e)\nabla^2\bar{u}_i + 2(\partial_j\nu_e)\bar{S}_{ij}$. Keep the full symmetric stress inside the divergence when $\nu_e$ varies in space; multiplying the Laplacian by the local viscosity, or taking the divergence of an unsymmetrized velocity gradient, omits part of the SGS force. compute_sgs_stress returns $\tau_{ij}=-2\nu_e\bar{S}_{ij}$; its contribution to momentum is $-\partial_j\tau_{ij}$.

Energy Cascade and Dissipation

SGS Dissipation

The rate at which energy is transferred from resolved to unresolved scales is:

\[\varepsilon_{sgs} = -\tau_{ij} \bar{S}_{ij} = 2\nu_e \bar{S}_{ij}\bar{S}_{ij} = \nu_e |\bar{S}|^2\]

where $|\bar{S}| = \sqrt{2\bar{S}_{ij}\bar{S}_{ij}}$ is the strain rate magnitude.

There is no extra factor of 2 in the final expression: the $\sqrt{2}$ in the definition of $|\bar{S}|$ already carries it. This is what sgs_dissipation returns — writing $2\nu_e|\bar{S}|^2$ would double-count the dissipation.

Available Models in Tarang.jl

Tarang.jl provides two eddy viscosity models:

Smagorinsky Model

Physical Basis

Joseph Smagorinsky (1963) proposed the first and most widely-used SGS model. It's based on dimensional analysis and mixing length theory.

Key assumptions:

1. The subgrid scales are in equilibrium (production = dissipation)
2. The characteristic length scale is proportional to the filter width $\Delta$
3. The characteristic velocity scale is $\Delta |\bar{S}|$

Mathematical Formulation

The eddy viscosity is computed as:

\[\nu_e = (C_s \Delta)^2 |\bar{S}|\]

where:

• \[C_s \approx 0.17\]

  is the Smagorinsky constant

• \[\Delta = (\Delta_x \Delta_y \Delta_z)^{1/3}\]

  is the effective filter width

• \[|\bar{S}| = \sqrt{2 \bar{S}_{ij} \bar{S}_{ij}}\]

  is the strain rate magnitude

Expanded form for 3D:

\[|\bar{S}| = \sqrt{2\left[\left(\frac{\partial \bar{u}}{\partial x}\right)^2 + \left(\frac{\partial \bar{v}}{\partial y}\right)^2 + \left(\frac{\partial \bar{w}}{\partial z}\right)^2\right] + \left(\frac{\partial \bar{u}}{\partial y} + \frac{\partial \bar{v}}{\partial x}\right)^2 + \left(\frac{\partial \bar{u}}{\partial z} + \frac{\partial \bar{w}}{\partial x}\right)^2 + \left(\frac{\partial \bar{v}}{\partial z} + \frac{\partial \bar{w}}{\partial y}\right)^2}\]

Usage in Tarang.jl

using Tarang

# Small runnable grid (use your simulation's actual grid size in production)
N = 16                      # Grid points per direction
L = 2π                      # Domain size
Δ = L / N                   # Grid spacing

sgs_model = SmagorinskyModel(
    C_s = 0.17,                    # Smagorinsky constant
    filter_width = (Δ, Δ, Δ),      # Grid spacing in each direction
    field_size = (N, N, N)         # Number of grid points
)

# Compute eddy viscosity from all 9 velocity-gradient components. These
# synthetic arrays make the snippet runnable; replace them with your gradients.
gradients = ntuple(i -> fill(i == 1 ? 1.0 : 0.0, N, N, N), 9)
compute_eddy_viscosity!(sgs_model, gradients...)

νₑ = get_eddy_viscosity(sgs_model)

Choosing the Smagorinsky Constant

The constant $C_s$ is not universal and depends on the flow:

Rule of thumb: Start with $C_s = 0.17$ and reduce if you observe over-dissipation.

Limitations

1. Over-dissipation: The model is always "on", even in laminar regions
2. No backscatter: Cannot represent energy transfer from small to large scales
3. Wall behavior: Needs damping functions near solid walls
4. Isotropic assumption: Assumes same behavior in all directions

Anisotropic Minimum Dissipation (AMD) Model

Motivation

The AMD model (Rozema et al., 2015) addresses key limitations of Smagorinsky:

1. Automatic switch-off: $\nu_e = 0$ in laminar/transitional regions
2. Anisotropic grids: Properly handles $\Delta_x \neq \Delta_y \neq \Delta_z$
3. Minimum dissipation: Adds only the dissipation needed for stability

Mathematical Formulation

The AMD model computes eddy viscosity as:

\[\nu_e = \max\left(0, -C \frac{\Delta_k^2 \frac{\partial \bar{u}_i}{\partial x_k} \frac{\partial \bar{u}_j}{\partial x_k} \bar{S}_{ij}}{\frac{\partial \bar{u}_m}{\partial x_n} \frac{\partial \bar{u}_m}{\partial x_n}}\right)\]

with summation over $i$, $j$ and $k$. The filter width $\Delta_k$ is indexed by the derivative direction $k$, not by the velocity-component index — that is what makes the model anisotropy-aware, and it is the term most easily got wrong. Tying $\Delta$ to the component index instead (e.g. $\hat{\delta}_{ij} = \Delta_i \delta_{ij}$) leaves only the diagonal strain contributing, and gives a different sign as well as a different magnitude.

Simplified form:

\[\nu_e = \max\left(0, -C \frac{\text{numerator}}{\text{denominator}}\right)\]

• Numerator: Measures alignment between velocity gradients and strain rate
• Denominator: Total velocity gradient magnitude squared
• max(0, ...): Ensures non-negative eddy viscosity

Physical Interpretation

The numerator can be negative (giving $\nu_e > 0$) when:

• Velocity gradients are aligned with the strain rate
• Energy is being transferred to smaller scales

The numerator is positive (giving $\nu_e = 0$) when:

• Flow is laminar or transitional
• No SGS dissipation is needed

Usage in Tarang.jl

sgs_model = AMDModel(
    C = 1/12,                      # Poincaré constant (for spectral methods)
    filter_width = (Δx, Δy, Δz),   # Can be different in each direction
    field_size = (Nx, Ny, Nz),
    clip_negative = true           # Ensure νₑ ≥ 0 (recommended)
)

compute_eddy_viscosity!(sgs_model,
    ∂u∂x, ∂u∂y, ∂u∂z,
    ∂v∂x, ∂v∂y, ∂v∂z,
    ∂w∂x, ∂w∂y, ∂w∂z
)

νₑ = get_eddy_viscosity(sgs_model)

Numerical Evaluation

The CPU and GPU implementations use the same scalar kernels. Before evaluating the contractions, AMD normalizes velocity gradients and, for diffusivity, scalar gradients independently. It restores the velocity scale afterward; the scalar-gradient scale cancels. This preserves the expected scaling across very small and large gradient magnitudes without an absolute epsilon cutoff. Zero gradients give zero closure coefficients, while nonfinite gradients remain visible as NaN. Clipping negative predictions still follows clip_negative.

Choosing the AMD Constant

The constant $C$ depends on the numerical discretization:

Scalar Transport (Buoyancy)

For buoyancy-driven flows, AMD also provides eddy diffusivity for scalar transport:

\[\kappa_e = \max\left(0, -C \frac{\Delta_k^2 \frac{\partial \bar{u}_i}{\partial x_k} \frac{\partial \bar{\theta}}{\partial x_k} \frac{\partial \bar{\theta}}{\partial x_i}}{\frac{\partial \bar{\theta}}{\partial x_n} \frac{\partial \bar{\theta}}{\partial x_n}}\right)\]

The inner sum runs over every velocity component $u_i$, so the method needs the complete velocity-gradient tensor (9 components in 3D, 4 in 2D) before the scalar gradients — not the gradient of a single velocity component:

# Compute eddy diffusivity for buoyancy/temperature
compute_eddy_diffusivity!(sgs_model,
    ∂u∂x, ∂u∂y, ∂u∂z,    # Gradients of u
    ∂v∂x, ∂v∂y, ∂v∂z,    # Gradients of v
    ∂w∂x, ∂w∂y, ∂w∂z,    # Gradients of w
    ∂b∂x, ∂b∂y, ∂b∂z     # Buoyancy/scalar gradients
)

κₑ = get_eddy_diffusivity(sgs_model)

Choosing Between Models

Summary Comparison

Complete Example: LES of Decaying Turbulence

This example builds a model, evaluates gradients, and refreshes νₑ for a Taylor–Green velocity field. Insert the momentum update in the indicated location.

using Tarang
using Statistics

# 1. Physical and Numerical Parameters

N = 16                      # Smoke-test resolution; use 128+ for production
L = 2π                      # Domain size [m]
Δ = L / N                   # Grid spacing [m]
ν = 1e-4                    # Molecular (kinematic) viscosity [m²/s]
dt = 0.001                  # Time step [s]
nsteps = 3                  # Smoke run; increase for a production simulation

# 2. Create the SGS Model

sgs = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (Δ, Δ, Δ),
    field_size = (N, N, N),
    architecture = CPU()       # GPU() keeps νₑ on the device
)

# 3. Setup Computational Domain

coords = CartesianCoordinates("x", "y", "z")
dist = Distributor(coords; mesh=(1, 1, 1))  # Single processor

xbasis = RealFourier(coords["x"]; size=N, bounds=(0.0, L))
ybasis = RealFourier(coords["y"]; size=N, bounds=(0.0, L))
zbasis = RealFourier(coords["z"]; size=N, bounds=(0.0, L))

domain = Domain(dist, (xbasis, ybasis, zbasis))

u = VectorField(dist, coords, "u", (xbasis, ybasis, zbasis))

# 4. Initialize with a Turbulent Velocity Field

# Taylor–Green initial condition
x, y, z = local_grids(dist, xbasis, ybasis, zbasis)
ensure_layout!(u, :g)
get_grid_data(u.components[1]) .=  sin.(x) .* cos.(y') .* cos.(reshape(z, 1, 1, :))
get_grid_data(u.components[2]) .= .-cos.(x) .* sin.(y') .* cos.(reshape(z, 1, 1, :))
get_grid_data(u.components[3]) .= 0.0
ensure_layout!(u, :c)

# 5. Helper Function: Compute All Velocity Gradients

function compute_velocity_gradients(u)
    # evaluate(grad(u)) returns a TensorField in component-major order:
    # (∂u∂x, ∂u∂y, ∂u∂z, ∂v∂x, ∂v∂y, ∂v∂z, ∂w∂x, ∂w∂y, ∂w∂z).
    G = evaluate(grad(u))
    return Tuple(get_grid_data(g) for g in G.components)
end

# 6. Per-Step SGS Update

for step in 1:nsteps
    grads = compute_velocity_gradients(u)

    compute_eddy_viscosity!(sgs, grads...)

    νₑ = get_eddy_viscosity(sgs)
    ν_eff = ν .+ νₑ  # Total viscosity = molecular + SGS, one value per grid point

    # Apply your momentum update here using the full symmetric stress:
    #   ∂ū/∂t + (ū·∇)ū = -∇p*/ρ + ∇·[2(ν + νₑ)S̄]
    #   S̄ = (∇ū + transpose(∇ū))/2
    #
    # Assemble the variable-viscosity term explicitly from ν_eff and S̄.

    if step == nsteps
        mean_νₑ = mean_eddy_viscosity(sgs)
        max_νₑ = max_eddy_viscosity(sgs)

        println("Step $step:")
        println("  Mean eddy viscosity: $(mean_νₑ)")
        println("  Max eddy viscosity:  $(max_νₑ)")
        println("  Ratio νₑ/ν (mean):   $(mean_νₑ/ν)")
    end
end

Complete Example: LES of Rayleigh-Bénard Convection

For buoyancy-driven turbulence we need both momentum and scalar closures. As above, only the SGS update is shown — the momentum and temperature updates are yours to write.

using Tarang

# Physical Parameters

Nx, Nz = 256, 128           # Grid resolution
Lx, Lz = 4.0, 1.0           # Domain size

Ra = 1e8                    # Rayleigh number
Pr = 1.0                    # Prandtl number
ν = sqrt(Pr / Ra)           # Kinematic viscosity
κ = ν / Pr                  # Thermal diffusivity
nsteps = 1000               # Number of time steps

# Domain and Fields

coords = CartesianCoordinates("x", "z")
dist   = Distributor(coords; dtype=Float64, device=CPU())
xbasis = RealFourier(coords["x"]; size=Nx, bounds=(0.0, Lx))  # periodic
zbasis = ChebyshevT(coords["z"];  size=Nz, bounds=(0.0, Lz))  # wall-bounded
domain = Domain(dist, (xbasis, zbasis))

u = Field(dist; name="u", bases=(xbasis, zbasis))   # horizontal velocity
w = Field(dist; name="w", bases=(xbasis, zbasis))   # vertical velocity
T = Field(dist; name="T", bases=(xbasis, zbasis))   # temperature

# (Initialise u, w and T with your own IC and boundary conditions.)

# Create AMD Model (recommended for RBC)

# One filter width per axis. `grid_spacing` reports the *smallest* spacing on
# each axis, which is the conservative choice on the Chebyshev-clustered z axis.
Δx, Δz = grid_spacing(domain)

sgs = AMDModel(
    C = 1/12,
    filter_width = (Δx, Δz),
    field_size = (Nx, Nz)
)

# Per-Step SGS Update

for step in 1:nsteps
    # Velocity gradients. grad() is lazy and not iterable; evaluate() realises
    # it, and for a scalar field the result has one component per axis.
    ∂u∂x, ∂u∂z = (get_grid_data(g) for g in evaluate(grad(u)).components)
    ∂w∂x, ∂w∂z = (get_grid_data(g) for g in evaluate(grad(w)).components)

    # Temperature gradients
    ∂T∂x, ∂T∂z = (get_grid_data(g) for g in evaluate(grad(T)).components)

    # --- Momentum closure ---
    # 2D needs the complete 2×2 velocity-gradient tensor, component-major.
    compute_eddy_viscosity!(sgs, ∂u∂x, ∂u∂z, ∂w∂x, ∂w∂z)
    νₑ = get_eddy_viscosity(sgs)
    ν_eff = ν .+ νₑ

    # --- Scalar (temperature) closure ---
    # The same velocity-gradient tensor first, then the scalar gradients.
    compute_eddy_diffusivity!(sgs, ∂u∂x, ∂u∂z, ∂w∂x, ∂w∂z, ∂T∂x, ∂T∂z)
    κₑ = get_eddy_diffusivity(sgs)
    κ_eff = κ .+ κₑ

    # ν_eff and κ_eff are plain arrays. Applying them to the momentum and
    # temperature equations is your code — see the scope note at the top.
end

Diagnostics and Analysis

SGS Energy Dissipation

The SGS dissipation rate tells you how much energy is being drained by the model:

# Get strain magnitude (computed during eddy viscosity calculation)
S_mag = sgs.strain_magnitude  # Only for Smagorinsky

# Or compute it yourself — note the broadcasting dots (sqrt is not vectorised)
S_mag = sqrt.(2 .* (S11.^2 .+ S22.^2 .+ S33.^2 .+ 2 .* S12.^2 .+ 2 .* S13.^2 .+ 2 .* S23.^2))

# SGS dissipation field: ε_sgs = νₑ |S̄|²  (no extra factor of 2 — |S̄| carries it)
ε_sgs = sgs_dissipation(sgs, S_mag)

# Domain-averaged SGS dissipation
ε_sgs_mean = mean_sgs_dissipation(sgs, S_mag)

Monitoring Model Behavior

# Check ratio of SGS to molecular viscosity
νₑ_mean = mean_eddy_viscosity(sgs)
println("νₑ/ν = $(νₑ_mean/ν)")

# If νₑ/ν >> 1: SGS model is dominant (typical for high Re)
# If νₑ/ν << 1: Flow is nearly DNS-resolved
# If νₑ/ν ~ 1: Well-resolved LES

# For AMD: check how often νₑ = 0
νₑ = get_eddy_viscosity(sgs)
fraction_zero = sum(νₑ .== 0) / length(νₑ)
println("Fraction with νₑ = 0: $(fraction_zero)")
# High fraction → flow is largely laminar/transitional

Tips for Successful LES

Check grid and timestep convergence, monitor SGS dissipation, and use the actual axis spacings for filter widths. See Troubleshooting for common symptoms; viscosity alone does not establish whether a run is resolved.

When LES May Not Be Appropriate

• Very low Reynolds numbers: DNS may be feasible
• Strongly anisotropic turbulence: May need special treatment
• Flows with strong backscatter: Standard models don't capture this
• Near-wall regions: May need wall models or finer grids

2D Flows

Both models support 2D simulations:

sgs_2d = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (Δx, Δy),
    field_size = (Nx, Ny)
)

# 2D velocity gradients (4 components)
compute_eddy_viscosity!(sgs_2d, ∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y)

Note on 2D turbulence: 2D turbulence has an inverse energy cascade (energy flows to large scales), which is fundamentally different from 3D. Standard SGS models may not be appropriate for 2D flows.

Choosing an Architecture (CPU / GPU)

Both models take an architecture keyword, defaulting to CPU(). It decides where the model's internal buffers (eddy_viscosity, strain_magnitude, and for AMD eddy_diffusivity) are allocated: plain Arrays on CPU(), CuArrays on GPU().

# CPU (the default)
sgs_cpu = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (Δ, Δ, Δ),
    field_size = (N, N, N),
    architecture = CPU()
)

# AMD takes the same keyword
amd_cpu = AMDModel(
    C = 1/12,
    filter_width = (Δx, Δy, Δz),
    field_size = (Nx, Ny, Nz),
    architecture = CPU()
)

GPU() requires CUDA.jl to be loaded first, and a CUDA device to be present — it throws otherwise:

using CUDA        # must be loaded before GPU() can be constructed
using Tarang

sgs_gpu = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (Δ, Δ, Δ),
    field_size = (N, N, N),
    architecture = GPU()
)

compute_eddy_viscosity! and compute_eddy_diffusivity! broadcast the same scalar kernels on CPU and GPU arrays. GPU models upload CPU gradients on each call; keep gradients on the device to avoid repeated transfers. CPU models reject GPU gradients. get_eddy_viscosity returns the model's array type (Array on CPU, CuArray on CUDA).

API Reference

See LES Models API for constructor keywords, gradient ordering, MPI array requirements, diagnostics, and utility functions.

Troubleshooting

References

Original Papers

1. Smagorinsky, J. (1963). "General circulation experiments with the primitive equations: I. The basic experiment." Monthly Weather Review, 91(3), 99-164.

2. Rozema, W., Bae, H. J., Moin, P., & Verstappen, R. (2015). "Minimum-dissipation models for large-eddy simulation." Physics of Fluids, 27(8), 085107. PDF

3. Abkar, M., Bae, H. J., & Moin, P. (2016). "Minimum-dissipation scalar transport model for large-eddy simulation of turbulent flows." Physical Review Fluids, 1(4), 041701. Link

Textbooks and Reviews

1. Pope, S. B. (2000). Turbulent Flows. Cambridge University Press. (Chapter 13: Large-Eddy Simulation)

2. Sagaut, P. (2006). Large Eddy Simulation for Incompressible Flows. Springer.

Related Documentation

• dedaLES AMD Documentation

See Also

• Solvers - Time integration methods
• Operators - Gradient and differential operators
• Stochastic Forcing - Forcing for turbulence
• API: LES Models - Complete API reference