LES Models API

The SGS models are array-level utilities: they consume grid-space velocity-gradient arrays and produce an eddy-viscosity array. They are not automatically coupled into an IVP — you evaluate the gradients, call the model, and apply the resulting stress yourself.

Usage

using Tarang

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

u = VectorField(domain, "u")
x, y = local_grids(dist, xbasis, ybasis)
ensure_layout!(u, :g)
get_grid_data(u.components[1]) .=  sin.(x) .* cos.(y')
get_grid_data(u.components[2]) .= .-cos.(x) .* sin.(y')
ensure_layout!(u, :c)

# grad(u) is a TensorField; its components are the velocity-gradient tensor in
# component-major order: (∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y) — exactly the order the models want.
G = evaluate(grad(u))
∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y = (get_grid_data(g) for g in G.components)

Δx, Δy = grid_spacing(domain)
model  = SmagorinskyModel(C_s=0.17, filter_width=(Δx, Δy), field_size=(16, 16))

compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y)

νₑ = get_eddy_viscosity(model)          # 16×16 Matrix{Float64}
mean_eddy_viscosity(model)              # 0.00352…
max_eddy_viscosity(model)               # 0.00891…

# deviatoric SGS stress τᵢⱼ = -2 νₑ S̄ᵢⱼ
S11 = ∂u∂x
S12 = 0.5 .* (∂u∂y .+ ∂v∂x)
S22 = ∂v∂y
τ11, τ12, τ22 = compute_sgs_stress(model, S11, S12, S22)

mean_sgs_dissipation(model, model.strain_magnitude)   # 0.00643…

Usage under MPI

The models hold plain per-rank arrays, so size the model to this rank's slab, not to the global grid. Under MPI get_grid_data returns a PencilArray whose size is already the local shape; parent gives the raw local storage:

G = evaluate(grad(u))
∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y = (parent(get_grid_data(g)) for g in G.components)

Δx, Δy = grid_spacing(domain)
model  = SmagorinskyModel(C_s=0.17, filter_width=(Δx, Δy),
                          field_size=size(∂u∂x))   # LOCAL slab: (16, 8) at np=2

compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y)

mean_eddy_viscosity(model)      # 0.00352… — global, Allreduced across ranks
max_eddy_viscosity(model)       # 0.00891… — global

mean_eddy_viscosity, max_eddy_viscosity and mean_sgs_dissipation each perform their own reduction over MPI.COMM_WORLD, so they return domain-global values: for the 16×16 case above they give the same numbers at np = 1, 2 and 4. Being collectives, they must be called on every rank — reaching for one inside an if rank == 0 block deadlocks. Pass global_reduce=false to skip the reduction and get this rank's slab value, which is safe to call from a subset of ranks:

# Deadlocks at np > 1 — only rank 0 enters the collective.
rank == 0 && @info "mean νₑ" mean_eddy_viscosity(model)

# Correct: every rank reduces, one rank prints.
ν̄ = mean_eddy_viscosity(model)
rank == 0 && @info "mean νₑ" ν̄

# Or opt out of the reduction entirely.
rank == 0 && @info "local νₑ" mean_eddy_viscosity(model; global_reduce=false)

get_eddy_viscosity(model) is, by contrast, this rank's slab — broadcast it against the other parent(...) arrays to build the stress locally.

Sizing the model to the global grid is a loud error rather than a silent wrong answer:

DimensionMismatch: Gradient array 1 has size (16, 8), expected (16, 16)

Types

Abstract Types

abstract type SGSModel end
abstract type EddyViscosityModel <: SGSModel end

SmagorinskyModel

SmagorinskyModel{T, N, A, Arch}

νₑ = (Cₛ Δ)² |S̄|

# Create model for 256³ domain with Δ = 2π/256
model = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (2π/256, 2π/256, 2π/256),
    field_size = (256, 256, 256)
)

# Create GPU model
model_gpu = SmagorinskyModel(
    C_s = 0.17,
    filter_width = (2π/256, 2π/256, 2π/256),
    field_size = (256, 256, 256),
    architecture = GPU()
)

# Compute eddy viscosity from velocity gradients
compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂u∂z, ∂v∂x, ∂v∂y, ∂v∂z, ∂w∂x, ∂w∂y, ∂w∂z)

# Access the result
νₑ = get_eddy_viscosity(model)

Classic Smagorinsky (1963) subgrid-scale model.

Type signature:

mutable struct SmagorinskyModel{T<:AbstractFloat, N,
                                A<:AbstractArray{T, N},
                                Arch<:AbstractArchitecture} <: EddyViscosityModel

The array parameter A is the type the model's internal buffers are stored in — an Array{T,N} on CPU(), a CuArray on GPU().

Fields:

Constructor:

SmagorinskyModel(;
    C_s = 0.17,
    filter_width::NTuple{N, Real},
    field_size::NTuple{N, Int},
    dtype = Float64,
    architecture = CPU()
)

AMDModel

AMDModel{T, N, A, Arch}

νₑ = max(0, νₑ†)

νₑ† = -C (Δₖ² ∂uᵢ/∂xₖ ∂uⱼ/∂xₖ Sᵢⱼ) / (∂uₘ/∂xₙ ∂uₘ/∂xₙ)

# Create AMD model for anisotropic grid
model = AMDModel(
    C = 1/12,  # Spectral method
    filter_width = (2π/256, 2π/256, 2π/64),  # Anisotropic
    field_size = (256, 256, 64)
)

# Create GPU AMD model
model_gpu = AMDModel(
    C = 1/12,
    filter_width = (2π/256, 2π/256, 2π/64),
    field_size = (256, 256, 64),
    architecture = GPU()
)

# Compute eddy viscosity
compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂u∂z, ∂v∂x, ∂v∂y, ∂v∂z, ∂w∂x, ∂w∂y, ∂w∂z)

Anisotropic Minimum Dissipation model (Rozema et al., 2015).

Type signature:

mutable struct AMDModel{T<:AbstractFloat, N,
                        A<:AbstractArray{T, N},
                        Arch<:AbstractArchitecture} <: EddyViscosityModel

Fields:

Constructor:

AMDModel(;
    C = 1/12,
    filter_width::NTuple{N, Real},
    field_size::NTuple{N, Int},
    clip_negative = true,
    dtype = Float64,
    architecture = CPU()
)

Eddy Viscosity Computation

computeeddyviscosity!

compute_eddy_viscosity!(model::AMDModel, velocity_gradients...)

compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y)

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

Compute eddy viscosity from velocity gradient components. The gradients are plain grid-space arrays of size model.field_size, passed in component-major order (all derivatives of u, then all of v, then all of w).

2D Signature:

compute_eddy_viscosity!(model, ∂u∂x, ∂u∂y, ∂v∂x, ∂v∂y)

3D Signature:

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

Returns: The eddy viscosity array model.eddy_viscosity

computeeddydiffusivity!

compute_eddy_diffusivity!(model::AMDModel, velocity_gradients..., scalar_gradients...)

κₑ = max(0, κₑ†),   κₑ† = -C · [ Σₖ δₖ² (∂ₖ uᵢ)(∂ₖ b)(∂ᵢ b) ] / [ (∂ₗ b)(∂ₗ b) ]

Compute eddy diffusivity for scalar transport (AMD model only). The AMD diffusivity (Abkar, Bae & Moin 2016, eq. 2.7) is a full double contraction over the scaled-gradient direction k and every velocity component i:

κₑ† = -C · [ Σₖ Δₖ² (∂ₖuᵢ)(∂ₖb)(∂ᵢb) ] / [ (∂ₗb)(∂ₗb) ]

so the method needs the complete velocity-gradient tensor (2D: 4 components, 3D: 9), not just the gradient of one velocity component. The scalar gradients follow.

2D Signature:

compute_eddy_diffusivity!(model::AMDModel,
    ∂u∂x, ∂u∂y,
    ∂v∂x, ∂v∂y,
    ∂b∂x, ∂b∂y
)

3D Signature:

compute_eddy_diffusivity!(model::AMDModel,
    ∂u∂x, ∂u∂y, ∂u∂z,
    ∂v∂x, ∂v∂y, ∂v∂z,
    ∂w∂x, ∂w∂y, ∂w∂z,
    ∂b∂x, ∂b∂y, ∂b∂z
)

Returns: The eddy diffusivity array model.eddy_diffusivity

Subgrid Stress

computesgsstress

compute_sgs_stress(model::EddyViscosityModel, strain_components...)

τᵢⱼᵈ = -2 νₑ S̄ᵢⱼ

Compute deviatoric SGS stress tensor τᵢⱼ = -2 νₑ S̄ᵢⱼ.

2D Signature:

compute_sgs_stress(model, S11, S12, S22) -> (τ11, τ12, τ22)

3D Signature:

compute_sgs_stress(model, S11, S12, S13, S22, S23, S33)
    -> (τ11, τ12, τ13, τ22, τ23, τ33)

Accessors

geteddyviscosity

get_eddy_viscosity(model::EddyViscosityModel)

Return the current eddy viscosity field (an Array on CPU(), a CuArray on GPU()).

get_eddy_viscosity(model::EddyViscosityModel) -> A <: AbstractArray{T, N}

geteddydiffusivity

get_eddy_diffusivity(model::AMDModel)

Return the current eddy diffusivity field (AMD only).

get_eddy_diffusivity(model::AMDModel) -> A <: AbstractArray{T, N}

getfilterwidth

get_filter_width(model::EddyViscosityModel)

Return the filter width tuple.

get_filter_width(model::EddyViscosityModel) -> NTuple{N, T}

Statistics

meaneddyviscosity

mean_eddy_viscosity(model::EddyViscosityModel; global_reduce=true)

Compute domain-averaged eddy viscosity.

mean_eddy_viscosity(model::EddyViscosityModel) -> T

maxeddyviscosity

max_eddy_viscosity(model::EddyViscosityModel; global_reduce=true)

Return maximum eddy viscosity in the domain.

max_eddy_viscosity(model::EddyViscosityModel) -> T

Dissipation Rate

sgs_dissipation

sgs_dissipation(model::EddyViscosityModel, strain_magnitude::AbstractArray)

εₛₛ = νₑ |S̄|²

Compute the SGS dissipation rate field: εₛₛ = νₑ |S̄|²

There is no extra factor of 2 here. The exact dissipation is εₛₛ = -τᵢⱼ S̄ᵢⱼ = 2 νₑ S̄ᵢⱼS̄ᵢⱼ, and the |S̄| = √(2 S̄ᵢⱼS̄ᵢⱼ) convention used by compute_eddy_viscosity! already carries it, so εₛₛ = νₑ |S̄|².

sgs_dissipation(model::EddyViscosityModel, strain_magnitude) -> A <: AbstractArray{T, N}

The natural argument is the strain magnitude the model just cached, model.strain_magnitude (Smagorinsky only — the AMD model does not cache one).

meansgsdissipation

mean_sgs_dissipation(model::EddyViscosityModel, strain_magnitude::AbstractArray)

Compute domain-averaged SGS dissipation rate.

mean_sgs_dissipation(model::EddyViscosityModel, strain_magnitude) -> T

Configuration

set_constant!

set_constant!(model::SmagorinskyModel, C_s::Real)

set_constant!(model::AMDModel, C::Real)

Update the model constant.

set_constant!(model::SmagorinskyModel, C_s::Real)
set_constant!(model::AMDModel, C::Real)

reset!

Reset eddy viscosity (and diffusivity for AMD) to zero.

reset!(model::EddyViscosityModel)
reset!(model::AMDModel)  # Also resets eddy_diffusivity

Exports

export SGSModel, EddyViscosityModel
export SmagorinskyModel, AMDModel
export compute_eddy_viscosity!, compute_eddy_diffusivity!
export compute_sgs_stress
export get_eddy_viscosity, get_eddy_diffusivity, get_filter_width
export mean_eddy_viscosity, max_eddy_viscosity
export reset!, set_constant!, set_filter_width!
export sgs_dissipation, mean_sgs_dissipation

Mutating a model

set_constant!(model, C) and set_filter_width!(model, Δ) validate their argument and keep the cached derived quantities (filter_width_sq, effective_delta) in step. Assigning model.filter_width directly is also honoured — the kernels re-derive Δ² and the geometric-mean Δ on every call — but it leaves those cached fields reading stale, so prefer the setter.

Index

• Tarang.AMDModel
• Tarang.SmagorinskyModel
• Tarang.compute_eddy_diffusivity!
• Tarang.compute_eddy_viscosity!
• Tarang.compute_sgs_stress
• Tarang.get_eddy_diffusivity
• Tarang.get_eddy_viscosity
• Tarang.get_filter_width
• Tarang.max_eddy_viscosity
• Tarang.mean_eddy_viscosity
• Tarang.mean_sgs_dissipation
• Tarang.set_constant!
• Tarang.sgs_dissipation