|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
Loading...
Searching...
No Matches
|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
Large-eddy simulation resolves the motions the grid can carry and models the rest. The modelling is concentrated in a single coefficient, and the whole of this page is about where that coefficient comes from, what the code actually computes, and which choices are yours.
The YAML keys are specified in Configuration Reference: Case YAML; this page explains the physics behind them. The runtime call order is in C Runtime Execution Map.
Filtering the incompressible Navier-Stokes equations at the grid scale leaves the nonlinear term with a residue that cannot be evaluated from the filtered field alone:
tau_ij = (u_i u_j)~ - u~_i u~_j
This is the subgrid stress: the momentum transfer done by eddies smaller than a cell. It involves the unfiltered velocity, which the simulation never had. Closing the equations means modelling it.
Smagorinsky's closure is the simplest one that works, and it is an eddy-viscosity ansatz: pretend the unresolved motions act like an extra viscosity that drains energy from the resolved field.
tau_ij ~= -2 (Cs Delta)^2 |S~| S~_ij = -2 C Delta^2 |S~| S~_ij
with Delta the filter width, S~_ij the resolved strain rate, and |S~| = sqrt(2 S~_ij S~_ij). The eddy viscosity is nu_t = C Delta^2 |S~|, added to the molecular viscosity in the momentum equations.
Everything then rests on the single number C. Note the notation carefully, because the code follows it: **C is the coefficient that multiplies Delta^2 |S~|, which is Cs^2 in the classical form, not Cs**. The stored field holds C. Reported diagnostics convert to Cs, which is the quantity with a value readers recognise.
constant_smagorinsky prescribes Cs and is done. It is cheap, predictable, and correct wherever a single number is appropriate.
It is not appropriate in most places. The right value differs between isotropic turbulence (Lilly's 0.16-0.17, valid only where the grid cutoff sits in an inertial range), a wall-bounded flow (where the coefficient must vanish at the wall and no constant does), a transitional region (where any nonzero value dissipates a flow that is not yet turbulent), and a laminar pocket (where the model should switch itself off and cannot).
The dynamic procedure exists to let the flow supply the number.
The insight (Germano et al. 1991) is that although the simulation cannot see below its grid, it can see between the grid scale and a coarser one. Take the already-filtered field and filter it again with a wider test filter of width Delta^ = r Delta, and three stresses come into view:
tau_ij — the stress at the grid filter. Unknown; this is what is being modelled.T_ij — the stress at the combined grid-and-test filter. Also unknown.L_ij = (u~_i u~_j)^ - u^_i u^_j — the Leonard stress, computable entirely from the resolved field.The identity connecting them is exact and purely algebraic, with no modelling in it:
L_ij = T_ij - (tau_ij)^
L_ij is the stress carried by the eddies between Delta and Delta^: the band the simulation can actually resolve and therefore measure. That is the whole trick.
Now apply the same Smagorinsky ansatz at both levels with the same coefficient:
tau_ij ~= -2 C Delta^2 |S~| S~_ij T_ij ~= -2 C Delta^^2 |S^| S^_ij
Substituting into the identity gives L_ij ~= C M_ij with
M_ij = -2 Delta^2 ( alpha |S^| S^_ij - (|S~| S~_ij)^ ), alpha = (Delta^/Delta)^2
M_ij are structurally different and must stay so. The first filters S~ and then forms the product, so it is a product of filtered quantities. The second forms the product |S~| S~_ij pointwise on the fine field and then filters the whole thing, so it is a filter of a product. Those are not equal, and their difference is the entire information content of the dynamic procedure. Written with both terms as the same expression, M_ij collapses to a fixed multiple of |S^| S^_ij, and the result is a well-scaled, dimensionally correct, flow-responsive coefficient that is not the Germano-Lilly coefficient. It cannot be cited as one. This is exactly the shape the implementation once had, which is why GermanoModelTensor is one function with one call site and why tests/c/test_les.c pins the two terms apart on a field where they provably differ.L_ij = C M_ij is six equations for one unknown. Lilly (1992) resolves the overdetermination by least squares, minimising the residual over the tensor components:
C = <L_ij M_ij> / <M_ij M_ij>
The angle brackets are not decoration. The derivation pulls C out of a filtering operation, which is only legitimate if C is constant over whatever set the brackets average. That assumption is the subject of the next section.
M_ij is made deviatoric before the contraction. The trace of L_ij is absorbed by the pressure in an incompressible flow and carries no information about the coefficient; projecting M_ij trace-free removes it from the contraction exactly, rather than relying on the discrete divergence being small enough to ignore.
5.2 Coefficient Averaging Entries selects the set the two contractions are averaged over, and it matters for three separate reasons.
The derivation requires it. Lilly's least squares assumes C is constant across the averaging set. A pointwise coefficient was already inside the test filter used to build M_ij, so the local form is internally inconsistent. Averaging over a region where C genuinely is uniform restores the consistency.
The denominator vanishes. Wherever the resolved strain is briefly small, M:M collapses and the pointwise quotient becomes a ratio of two small noisy numbers. Averaging over a set large enough to contain structure fixes this at the source rather than guarding it after the fact.
Backscatter should cancel where it physically does. Negative L:M is real. Averaged over a homogeneous set, forward and reverse transfer combine the way the physics intends. Clipped pointwise, only the negative tail is removed, and the mean dissipation comes out biased upward.
Two implementation properties are worth knowing. The sums are volume-weighted, so a refined region does not dominate an average merely by holding more cells. And the numerator and denominator are averaged separately before dividing: the mean of the quotients is a different number from the quotient of the means, and it is the second that the closure defines. tests/c/test_les.c separates them explicitly.
averaging.mode: homogeneous takes its directions from the case's periodicity unless they are named explicitly. This is not a shortcut; periodicity is already the declaration that the flow repeats along an axis, which is the same claim homogeneity makes.
The closure reads the resolved per-axis flags rather than re-inspecting the boundary faces. The two cannot disagree: the loader derives the flags from the boundary files and rejects any case whose opposite faces, or whose blocks, do not agree on periodicity. Those flags are also what the DMDA itself was built from and what the shared spatial target consults, so there is one answer to the question rather than three places that each work it out. It means the two canonical cases need no extra configuration:
Cs(t) is then a single readable curve rather than a noisy field.Cs(y) profile. The coefficient falls toward the wall on its own, which is the behaviour a constant coefficient cannot reproduce.Naming a direction the case does not declare periodic is accepted with a warning, since homogeneity is a claim only the user can make.
A negative C means nu_t < 0: energy flowing from the unresolved scales back to the resolved ones. Backscatter is physical, and in a locally averaged dynamic model it occurs at a large fraction of points in developed turbulence.
It is also destabilising if unbounded, so something has to limit it. The three modes at 5.3 Coefficient Limiting Entries differ in what they limit and at what cost.
The important design point is where the bound is applied. The clipping modes bound the coefficient's sign, which removes backscatter entirely and biases the mean dissipation. clipping.mode: none instead keeps the signed coefficient and bounds the total viscosity:
nu + nu_t >= min_viscosity_ratio * nu
That constrains the quantity which actually has to stay positive for the momentum operator to remain well posed, and it leaves the physics alone everywhere the bound is not active. It is the better guard where backscatter matters, and it pairs naturally with homogeneous averaging, where the averaged coefficient goes negative only when the whole homogeneous set is backscattering on balance.
The Smagorinsky closures have no wall treatment of their own. The coefficient is whatever the dynamic procedure measures or the constant model prescribes, right down to the first cell, and nothing damps it toward the wall. wale falls toward a wall by construction; see 5.2 Coefficient-Free Models: WALE and Vreman.
A wall function is configured separately at turbulence.wall_function, and it does reach this closure. Contra2Cart rebuilds the Cartesian field from the contravariant one at the top of each timestep, discarding the previous solve's near-wall correction, so the wall model is re-applied before the strain rates are formed. Without that the eddy viscosity would be built from the uncorrected near-wall field while the momentum equation used the corrected one, and on the coarse grids that motivate a wall model the uncorrected strain is under-predicted.
The correction is an imposed analytic profile, so the Germano identity evaluated at the cell it overwrites would measure that profile rather than resolved turbulence. Those cells are therefore excluded from both contraction sums whenever a wall model is active. The exclusion is not a preference: the data there is not a measurement, and averaging it in biases whatever set it belongs to. Which set depends on the mode - the whole domain under global, the wall plane under homogeneous retaining the wall-normal direction, and that cell alone under local - so no mode escapes it. A cell left with no data takes a zero coefficient, which is the right answer where the wall model is already supplying the stress.
averaging.mode: homogeneous does give a channel a wall-normal coefficient profile, so the coefficient can fall toward the wall on its own where the resolved field supports it. That is not a substitute for a wall treatment on a grid too coarse to resolve the sublayer.
Both models build the eddy viscosity from invariants of the resolved velocity gradient g_ij = du_i/dx_j, chosen so that it vanishes where there is nothing to model. Neither needs a test filter, an averaging set, clipping, or a coefficient field, which is what makes them usable on geometry with no homogeneous direction, where the dynamic procedure falls back to a noisy pointwise coefficient.
WALE (Nicoud & Ducros 1999) uses the traceless symmetric part of the squared gradient, S^d_ij = (g_ik g_kj + g_jk g_ki)/2 - delta_ij g_kk^2 / 3:
S^d is zero in pure shear, so a laminar shear layer gets no eddy viscosity, and near a wall it scales so that nu_t falls as the cube of wall distance, the behaviour a Smagorinsky model needs a damping function to imitate. It uses one scalar width Delta.
Vreman (2004) uses alpha_ij = du_j/dx_i and a tensor that weights each direction by its own spacing:
Delta_m alpha_mi is the change in u_i across the cell in direction m. On a curvilinear grid PICurv measures that change along the cell's own edge, edge_m . grad(u_i), with the edges from ComputeCellEdgeVectors, so the model keeps Vreman's per-direction resolution on a stretched, turned cell and is unchanged by the grid's orientation. B_beta vanishes for any locally two-component gradient, including pure shear. This per-direction weighting is the property that matters on cells much longer in one direction than the others: every scalar filter_width, corrected or not, must describe such a cell with a single length.
Both are purely dissipative and, because they respond to resolved gradients, small in a laminar or transitional stretch. That is correct physics, but it means neither supplies the damping a non-dissipative convection scheme may need before resolved turbulence develops.
Two widths enter, and both are configurable because neither has a universally right value on a curvilinear grid.
The grid filter width Delta is derived per cell by 5.1 Grid Filter Width Entries. The default cube root of the cell volume is exact for a cube and progressively optimistic as a cell is stretched, because a geometric mean is dominated by the short directions. max_edge errs the other way: on a wall-resolved mesh a wall-adjacent cell takes the streamwise spacing as its width, which inflates its eddy viscosity by the square of the aspect ratio. It suits wall-modelled or hybrid treatments, whose near-wall cells are meant to be coarse, rather than a mesh built to resolve the wall layer. Both non-default widths are measured along the cell's own grid directions, so an identical cell gets the same width in a straight section and in a bend.
The test-to-grid ratio Delta^/Delta sets alpha and defaults to 2.0, giving alpha = 4. The configured value is the ratio in each direction the kernel filters, and Delta is the cube root of the three directional widths. The box filters all three, so alpha is the ratio squared. simpson_ik filters only xi and zeta and leaves eta at the grid width, so alpha = ratio^(4/3), about 2.52 at the default. Squaring the ratio there, as the procedure once did, halved the dynamic coefficient on the HOM02 benchmark. The value is exposed because the effective width of the discrete three-point stencil is not exactly twice the grid spacing: matching second moments for a top-hat gives sqrt(6) ~= 2.449. Varying it is a standard sensitivity study, and the ingress rejects any value at or below 1, where the two filters coincide and the procedure has nothing to measure.
gradient_model.enabled adds the Clark tensor-diffusivity term to the viscous flux. It is the leading term of a Taylor expansion of the Leonard stress, and it is structurally different from an eddy viscosity: it is not purely dissipative and it does not require a coefficient to be measured.
The stress is tau_ij = sum_k (Delta_k^2 / 12) du_i/dx_k du_j/dx_k, with Delta_k the cell's spacing along grid direction k. The implementation forms it from the velocity differences between neighbouring cells along each grid direction, which already carry that spacing, so it needs no separate width and uses none: filter_width does not affect this term. tests/c/test_poisson_rhs.c case clark-gradient-model-scales-with-filter-width-squared pins both the power of Delta and the 1/12 coefficient, on cells whose spacing is not 1, where an error in that power would otherwise cancel.
Combined with a Smagorinsky closure this is the classical mixed model — the gradient term supplies the correct stress structure while the eddy viscosity supplies the dissipation the gradient term lacks. The two are configured independently because they are independent terms.
With diagnostics.enabled, the solver appends one row per step (or per cadence steps) to <run.analysis.metrics>/les_coefficient.csv. Three groups of columns answer different questions.
Is the coefficient right? cs_effective is the whole-domain sqrt(<L:M>/<M:M>), computed the same way regardless of averaging mode so it stays comparable across modes. Lilly's analysis puts it near 0.16-0.17 for decaying isotropic turbulence with the cutoff in an inertial range; on the HOM02 benchmark at 64^3 with the box test filter it settles at 0.193-0.196 under homogeneous averaging and 0.181-0.183 under local averaging (10. Status and Evidence), so read a value in that range as the expected one on comparable grids. cs_mean, coefficient_rms, coefficient_min, and coefficient_max describe the spread; under global averaging the spread should be zero by construction.
How much is the model doing? nu_t_mean, nu_t_max, and nu_t_over_nu_mean give the eddy-viscosity level. k_sgs_mean is the modelled subgrid kinetic energy from the Yoshizawa relation k_sgs = 2 C_I Delta^2 |S~|^2. Compared against the resolved turbulent kinetic energy from Field Statistics it measures how much of the energy the grid resolves — Pope's criterion, which asks for the resolved share to exceed roughly 80%. The two halves come from different pipelines on purpose: resolved turbulent energy needs the mean subtracted correctly, which the statistics windows do and the LES module has no business doing.
What is the limiting costing? backscatter_fraction and limited_fraction are the volume fractions that were negative, and that the clip modified, before clipping. No stored field preserves that state, which is why these columns exist. A max_cs that never binds reads near zero.
Nu_t and CS come from Field Statistics instead, which can target them like any other catalogued field.The module is built from kernels rather than from one procedure, so a new closure reuses the parts it needs instead of restating them. The public ones are declared in include/les.h:
include/statistics_target.h, not here, because averaging fields over a spatial domain is not an LES concern; it turns a local coefficient into a plane- or volume-averaged one and needs no separate code path per mode;Another eddy-viscosity closure (WALE, Vreman, the sigma model) is one function consuming the strain rate and the filter width and producing C, plus the selector wiring in Modular Selector Extension Guide. Dispatch is an enum and a branch in FlowSolver, the same pattern the momentum solvers use; there is no registry to register with.
A different averaging strategy, such as the Lagrangian dynamic model of Meneveau, Lund and Cabot (1996), replaces PicurvSpatialRatioAverage with a relaxation transported along pathlines. It is not implemented and no placeholder exists for it.
A transported-PDF closure such as a velocity filtered density function is a different kind of model: it evolves particles rather than adding momentum diffusion, and the Lagrangian machinery it needs is already in the tree. What it needs from here is nu_t, Delta, and the subgrid kinetic energy — the first two are fields and the third is SubgridKineticEnergy. No accessor is provided in advance, because an interface with no caller is a guess.
All four models are supported, together with the homogeneous, global and local averaging modes, the clamp, clip_negative and none clipping modes, the volume_weighted_box test filter, and all four widths. geometric_mean rests on hom02-geometric-mean-equivalence-2026-10-09, which shows it reproduces cube_root_volume on orthogonal cells; it has not been run on skewed cells, where it differs. simpson_ik is supported for planes homogeneous in xi and zeta; the dynamic procedure weights it as a two-direction filter (section 6), after the squared ratio it once used halved the coefficient. The unit coverage: tests/c/test_les.c checks the model tensor against its closed form on constant strain, pins the filtered product apart from the product of filtered factors, verifies that the procedure returns exactly zero on uniform flow, checks that global averaging gives one coefficient per block, and requires the Simpson coefficient to scale as 1/(ratio^(4/3) - 1); tests/c/test_mpi_kernels.c checks that the averaging reduction does not depend on the decomposition.
Five measurements add to it.
hom02-les-execution-2026-10-09): seventeen runs from Wray's 512^3 DNS field at 64^3, 32^3 and 64 x 64 x 32. Every model removes the energy pile-up that the unmodelled run builds at the grid cutoff (mean log10 error +0.72), to between +0.01 (dynamic, homogeneous) and +0.14 (Vreman). The dynamic coefficient settles without drift at 0.193-0.196 with homogeneous averaging and 0.181-0.183 with local averaging - above the 0.16-0.17 this section once anticipated, which assumed a spectral cutoff. global equals homogeneous to every printed digit in the periodic box, and on 2:1 cells the eddy viscosity orders with the width as the definitions require.hom02-les-accuracy-2026-10-09, not met): against criteria fixed before the runs, no model keeps the resolved energy decay within 10% of the DNS at all times. Every model runs 18-29% high in the first 1.5 time units, which coincides with starting from a truncated DNS field at its peak-dissipation instant with no subgrid energy; by the end the constant, Vreman and WALE models are 13-16% low, homogeneous dynamic 10% low, and local dynamic within 2%. Low-band spectra are within 0.1 RMS in log10 for every model except Vreman (0.107).wmles-channel-retau1000-werner-2026-10-09): dynamic Smagorinsky with homogeneous averaging and the Werner wall model reproduces Lee & Moser (2015) within criteria fixed beforehand - bulk-unit mean velocity to 1.1% RMS, outer-layer resolved Reynolds shear stress to 0.91-1.02 - with Cs 0.133.hom02-simpson-ik-2026-10-09, wmles-channel-retau1000-simpson-2026-10-09): on HOM02 its former weighting gave Cs 0.091 against the box filter's 0.194 and a +0.39 cutoff pile-up; corrected, 0.168 and +0.12. On the channel the corrected filter meets every criterion, with outer u' 6% above Lee & Moser against 16% before, and 20% less eddy viscosity than the box filter.wmles-channel-retau1000-werner-window2-2026-10-09): the same flow over t = 2200-2800 reproduces every profile criterion, with u_tau 3.2% high, just outside the 3% bound.humphrey-turbulent-bend-wmles-2026-10-09): global averaging through a 90-degree bend at Re_D = 40,000, 17,000 steps without a failed step; its comparison with measured profiles is exploratory.What would change this: a refined-grid rerun of the channel on the current build, a wall-resolved reference case, and a run on skewed cells, where geometric_mean differs from cube_root_volume.
The eddy viscosity Nu_t and the dynamic coefficient CS are written with every checkpoint. The constant, Vreman and WALE models keep no other state and recompute nu_t from the restored velocity on the first step. The dynamic procedure recomputes its coefficient on its own cadence, so a restart continues from the checkpointed coefficient until the next update rather than from zero. The wall-modelled channel above was a restart of this kind, from step 50000 of an earlier run, with the statistics window reset.
1.9.8