les.h File Reference

#include "variables.h"
#include "logging.h"
#include "Metric.h"
#include "setup.h"
#include "Filter.h"
#include "Boundaries.h"

Include dependency graph for les.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
SymTensor	SymTensorSelfOuter (Cmpnts v)
	Forms a symmetric tensor from a vector's outer product with itself.
SymTensor	SymTensorCombine (PetscReal a, SymTensor x, PetscReal b, SymTensor y)
	Forms the linear combination a*x + b*y.
PetscReal	SymTensorTrace (SymTensor t)
	Returns the trace t_kk.
SymTensor	SymTensorDeviator (SymTensor t)
	Removes the isotropic part, returning t_ij - (1/3) delta_ij t_kk.
PetscReal	SymTensorContract (SymTensor a, SymTensor b)
	Contracts two symmetric tensors as a_ij b_ij.
Cmpnts	SymTensorTimesVector (SymTensor t, Cmpnts v)
	Applies a symmetric tensor to a vector, returning t_ij v_j.
PetscReal	SymTensorNormSq (SymTensor t)
	Returns the squared Frobenius norm t_ij t_ij.
PetscErrorCode	StrainRateFromGradients (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, SymTensor *strain, PetscReal *magnitude)
	Builds the strain-rate tensor and its magnitude from a velocity gradient.
PetscErrorCode	ComputeCellFilterWidth (LESFilterWidthModel model, PetscReal aj, Cmpnts csi, Cmpnts eta, Cmpnts zet, PetscReal *delta)
	Computes one cell's grid filter width under the selected width model.
SymTensor	LeonardStress (Cmpnts velocity_filtered, SymTensor velocity_product_filtered)
	Forms the Leonard stress L_ij = (u_i u_j)^ - u^_i u^_j.
SymTensor	GermanoModelTensor (PetscReal delta, PetscReal alpha, PetscReal strain_magnitude_filtered, SymTensor strain_filtered, SymTensor strain_product_filtered)
	Forms the deviatoric Germano model tensor M_ij.
PetscErrorCode	ResolveLESAveragingDirections (UserCtx *user, PetscBool direction[3])
	Resolves which logical directions the dynamic coefficient is averaged over.
PetscReal	ClipModelCoefficient (PetscReal coefficient, const LESConfig *config, PetscBool *limited)
	Applies the configured admissible range to one model coefficient.
PetscReal	EddyViscosityFromCoefficient (PetscReal coefficient, PetscReal delta, PetscReal strain_magnitude, PetscReal molecular_viscosity, PetscReal min_viscosity_ratio)
	Builds the eddy viscosity from a model coefficient.
PetscReal	WALEEddyViscosity (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, PetscReal delta, PetscReal coefficient)
	WALE eddy viscosity (Nicoud & Ducros 1999) from the resolved velocity gradient.
PetscReal	VremanEddyViscosity (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, const Cmpnts edges[3], PetscReal coefficient)
	Vreman eddy viscosity (Vreman 2004), weighting each grid direction by its own spacing.
PetscReal	SubgridKineticEnergy (PetscReal yoshizawa_ci, PetscReal delta, PetscReal strain_magnitude)
	Returns the modelled subgrid kinetic energy at a cell.
PetscErrorCode	ComputeSmagorinskyConstant (UserCtx *user)
	Computes the dynamic Smagorinsky coefficient field for one block.
PetscErrorCode	ComputeEddyViscosityLES (UserCtx *user)
	Computes the turbulent eddy viscosity for one block.
PetscErrorCode	LogLESDiagnostics (UserCtx *user)
	Appends one row of LES coefficient statistics to the run's log directory.

Function Documentation

SymTensorSelfOuter()

SymTensor SymTensorSelfOuter

(

Cmpnts

v

)

Forms a symmetric tensor from a vector's outer product with itself.

Used to build the resolved velocity products the Leonard stress subtracts.

Parameters

[in]

v

Vector to square.

Returns

The tensor v_i v_j.

Forms a symmetric tensor from a vector's outer product with itself.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorSelfOuter()

Definition at line 144 of file les.c.

{
    SymTensor t;
 
    t.xx = v.x * v.x;
    t.xy = v.x * v.y;
    t.xz = v.x * v.z;
    t.yy = v.y * v.y;
    t.yz = v.y * v.z;
    t.zz = v.z * v.z;
    return t;
}

Here is the caller graph for this function:

SymTensorCombine()

SymTensor SymTensorCombine	(	PetscReal	a,
		SymTensor	x,
		PetscReal	b,
		SymTensor	y
	)

Forms the linear combination a*x + b*y.

One primitive covers addition, subtraction, and scaling, so tensor algebra at the call sites reads as algebra instead of six repeated component expressions.

Parameters

[in]	a	Scalar multiplying x.
[in]	x	First tensor.
[in]	b	Scalar multiplying y.
[in]	y	Second tensor.

Returns

The combined tensor.

Forms the linear combination a*x + b*y.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorCombine()

Definition at line 165 of file les.c.

{
    SymTensor t;
 
    t.xx = a * x.xx + b * y.xx;
    t.xy = a * x.xy + b * y.xy;
    t.xz = a * x.xz + b * y.xz;
    t.yy = a * x.yy + b * y.yy;
    t.yz = a * x.yz + b * y.yz;
    t.zz = a * x.zz + b * y.zz;
    return t;
}

Here is the caller graph for this function:

SymTensorTrace()

PetscReal SymTensorTrace

(

SymTensor

t

)

Returns the trace t_kk.

Parameters

[in]

t

Tensor to trace.

Returns

Sum of the three diagonal components.

Returns the trace t_kk.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorTrace()

Definition at line 186 of file les.c.

{
    return t.xx + t.yy + t.zz;
}

Here is the caller graph for this function:

SymTensorDeviator()

SymTensor SymTensorDeviator

(

SymTensor

t

)

Removes the isotropic part, returning t_ij - (1/3) delta_ij t_kk.

The Germano model tensor must be deviatoric so that the trace of the Leonard stress, which the incompressible pressure absorbs, cannot enter the contraction through discrete divergence error.

Parameters

[in]

t

Tensor to project.

Returns

The trace-free part of t.

Removes the isotropic part, returning t_ij - (1/3) delta_ij t_kk.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorDeviator()

Definition at line 199 of file les.c.

{
    const PetscReal isotropic = SymTensorTrace(t) / 3.0;
 
    t.xx -= isotropic;
    t.yy -= isotropic;
    t.zz -= isotropic;
    return t;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SymTensorContract()

PetscReal SymTensorContract	(	SymTensor	a,
		SymTensor	b
	)

Contracts two symmetric tensors as a_ij b_ij.

Only six components are stored, so the three off-diagonal products are counted twice. Centralising that factor is the point: written out at each call site it is a silent error whenever one term is missed.

Parameters

[in]	a	First tensor.
[in]	b	Second tensor.

Returns

The full nine-term contraction.

Contracts two symmetric tensors as a_ij b_ij.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorContract()

Definition at line 217 of file les.c.

{
    // The three off-diagonal components each stand for two entries of the full tensor.
    return a.xx * b.xx + a.yy * b.yy + a.zz * b.zz +
           2.0 * (a.xy * b.xy + a.xz * b.xz + a.yz * b.yz);
}

Here is the caller graph for this function:

SymTensorTimesVector()

Cmpnts SymTensorTimesVector	(	SymTensor	t,
		Cmpnts	v
	)

Applies a symmetric tensor to a vector, returning t_ij v_j.

The three stored off-diagonal components each stand for two entries of the full tensor, so a caller writing this out by hand has six chances to reference the wrong one. Centralising it is the same reason SymTensorContract() exists.

Parameters

[in]	t	Tensor to apply.
[in]	v	Vector to apply it to.

Returns

The resulting vector.

Applies a symmetric tensor to a vector, returning t_ij v_j.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorTimesVector()

Definition at line 232 of file les.c.

{
    Cmpnts result;
 
    result.x = t.xx * v.x + t.xy * v.y + t.xz * v.z;
    result.y = t.xy * v.x + t.yy * v.y + t.yz * v.z;
    result.z = t.xz * v.x + t.yz * v.y + t.zz * v.z;
    return result;
}

Here is the caller graph for this function:

SymTensorNormSq()

PetscReal SymTensorNormSq

(

SymTensor

t

)

Returns the squared Frobenius norm t_ij t_ij.

Parameters

[in]

t

Tensor to measure.

Returns

The self-contraction of t.

Returns the squared Frobenius norm t_ij t_ij.

Full API contract is documented with the header declaration in include/les.h.

See also

SymTensorNormSq()

Definition at line 250 of file les.c.

{
    return SymTensorContract(t, t);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StrainRateFromGradients()

PetscErrorCode StrainRateFromGradients	(	Cmpnts	dudx,
		Cmpnts	dvdx,
		Cmpnts	dwdx,
		SymTensor *	strain,
		PetscReal *	magnitude
	)

Builds the strain-rate tensor and its magnitude from a velocity gradient.

S_ij = 0.5 (du_i/dx_j + du_j/dx_i) and |S| = sqrt(2 S_ij S_ij). This is the single definition of both quantities; the dynamic procedure and the eddy-viscosity assembly consume it rather than each expanding the algebra themselves.

Parameters

[in]	dudx	Gradient of the x velocity component.
[in]	dvdx	Gradient of the y velocity component.
[in]	dwdx	Gradient of the z velocity component.
[out]	strain	Resolved strain-rate tensor; may be NULL if only the magnitude is wanted.
[out]	magnitude	Strain-rate magnitude |S|; may be NULL.

Returns

PetscErrorCode 0 on success.

Builds the strain-rate tensor and its magnitude from a velocity gradient.

Full API contract is documented with the header declaration in include/les.h.

See also

StrainRateFromGradients()

Definition at line 267 of file les.c.

{
    SymTensor s;
 
    PetscFunctionBeginUser;
 
    s.xx = dudx.x;
    s.xy = 0.5 * (dudx.y + dvdx.x);
    s.xz = 0.5 * (dudx.z + dwdx.x);
    s.yy = dvdx.y;
    s.yz = 0.5 * (dvdx.z + dwdx.y);
    s.zz = dwdx.z;
 
    if (strain) *strain = s;
    // |S| = sqrt(2 S_ij S_ij), with the off-diagonal doubling handled by the contraction.
    if (magnitude) *magnitude = PetscSqrtReal(2.0 * SymTensorNormSq(s));
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeCellFilterWidth()

PetscErrorCode ComputeCellFilterWidth	(	LESFilterWidthModel	model,
		PetscReal	aj,
		Cmpnts	csi,
		Cmpnts	eta,
		Cmpnts	zet,
		PetscReal *	delta
	)

Computes one cell's grid filter width under the selected width model.

LES_FILTER_WIDTH_CUBE_ROOT_VOLUME uses the cell volume alone and so ignores anisotropy; the other models take the cell's extents along its own grid directions from ComputeCellDirectionalExtents, so they follow the cell's shape and are unchanged by its orientation.

Parameters

[in]	model	Width model to apply.
[in]	aj	Inverse cell volume (the Jacobian determinant) at the cell.
[in]	csi	Xi-direction metric vector at the cell.
[in]	eta	Eta-direction metric vector at the cell.
[in]	zet	Zeta-direction metric vector at the cell.
[out]	delta	Resolved filter width.

Returns

PetscErrorCode 0 on success, or PETSC_ERR_ARG_OUTOFRANGE for a non-positive Jacobian.

Computes one cell's grid filter width under the selected width model.

Full API contract is documented with the header declaration in include/les.h.

See also

ComputeCellFilterWidth()

Definition at line 296 of file les.c.

{
    double dx = 0.0, dy = 0.0, dz = 0.0;
 
    PetscFunctionBeginUser;
    PetscCheck(delta != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Filter width destination cannot be NULL.");
    PetscCheck(aj > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Cell Jacobian must be positive to define a filter width; received %g.", (double)aj);
 
    if (model == LES_FILTER_WIDTH_CUBE_ROOT_VOLUME) {
        // Cell volume alone: exact for a cube, and increasingly optimistic as the
        // cell is stretched, because it averages the three extents geometrically.
        *delta = PetscPowReal(1.0 / aj, 1.0 / 3.0);
        PetscFunctionReturn(0);
    }
 
    // The cell's own extents, not the Cartesian components of its diagonal: a filter
    // width that changed when the same cell was rotated would hand a bend nine times the
    // eddy viscosity of an identical cell in a straight section.
    PetscCall(ComputeCellDirectionalExtents(aj, csi, eta, zet, &dx, &dy, &dz));
 
    switch (model) {
    case LES_FILTER_WIDTH_MAX_EDGE:
        // The largest unresolved scale the cell can carry; conservative on stretched grids.
        *delta = PetscMax(dx, PetscMax(dy, dz));
        break;
    case LES_FILTER_WIDTH_SCOTTI: {
        // Scotti, Meneveau & Lilly (1993): the cube root of the extents, multiplied by
        // f(a1, a2), where a1 and a2 are the two shorter extents over the longest. f is 1
        // for a cube and grows with the cell's distortion. It assumes the cut-off lies in
        // the inertial range in every direction, which a very long cell violates.
        const double longest = PetscMax(dx, PetscMax(dy, dz));
        const double shortest = PetscMin(dx, PetscMin(dy, dz));
        const double middle = dx + dy + dz - longest - shortest;
        const double ln_a1 = log(shortest / longest), ln_a2 = log(middle / longest);
        const double shape = cosh(sqrt((4.0 / 27.0) *
                                       (ln_a1 * ln_a1 - ln_a1 * ln_a2 + ln_a2 * ln_a2)));
        *delta = PetscPowReal(PetscMax(dx * dy * dz, LES_EPSILON), 1.0 / 3.0) * shape;
        break;
    }
    case LES_FILTER_WIDTH_GEOMETRIC_MEAN:
    default:
        *delta = PetscPowReal(PetscMax(dx * dy * dz, LES_EPSILON), 1.0 / 3.0);
        break;
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LeonardStress()

SymTensor LeonardStress	(	Cmpnts	velocity_filtered,
		SymTensor	velocity_product_filtered
	)

Forms the Leonard stress L_ij = (u_i u_j)^ - u^_i u^_j.

L_ij is the stress carried by the scales between the grid filter and the test filter. It is computable entirely from the resolved field, which is what makes the dynamic procedure possible.

Parameters

[in]	velocity_filtered	Test-filtered resolved velocity u^_i.
[in]	velocity_product_filtered	Test filter of the products (u_i u_j)^.

Returns

The Leonard stress tensor.

Forms the Leonard stress L_ij = (u_i u_j)^ - u^_i u^_j.

Full API contract is documented with the header declaration in include/les.h.

See also

LeonardStress()

Definition at line 359 of file les.c.

{
    // L_ij = (u_i u_j)^ - u^_i u^_j : the stress carried between the two filter widths.
    return SymTensorCombine(1.0, velocity_product_filtered,
                            -1.0, SymTensorSelfOuter(velocity_filtered));
}

Here is the call graph for this function:

Here is the caller graph for this function:

GermanoModelTensor()

SymTensor GermanoModelTensor	(	PetscReal	delta,
		PetscReal	alpha,
		PetscReal	strain_magnitude_filtered,
		SymTensor	strain_filtered,
		SymTensor	strain_product_filtered
	)

Forms the deviatoric Germano model tensor M_ij.

M_ij = -2 Delta^2 ( alpha |S^| S^_ij - (|S| S_ij)^ ), made trace-free.

The two terms are structurally different and must stay so. The first is the model evaluated on the test-filtered strain, the product of separately filtered quantities. The second is the grid-level model formed pointwise and then test-filtered, that is, the filter of a product. Their difference is the entire information content of the dynamic procedure: collapsing them into the same expression yields a well-scaled, flow-responsive, and wrong coefficient that is not the Germano-Lilly coefficient and cannot be cited as one.

Parameters

[in]	delta	Grid filter width at the cell.
[in]	alpha	Squared test-to-grid filter width ratio.
[in]	strain_magnitude_filtered	Test-filtered strain magnitude |S|^.
[in]	strain_filtered	Test-filtered strain tensor S^_ij.
[in]	strain_product_filtered	Test filter of the product (|S| S_ij)^.

Returns

The deviatoric model tensor.

Forms the deviatoric Germano model tensor M_ij.

Full API contract is documented with the header declaration in include/les.h.

See also

GermanoModelTensor()

Definition at line 374 of file les.c.

{
    // Term one: the model evaluated on the test-filtered strain. Both factors were
    // filtered separately, then multiplied.
    const SymTensor test_scale = SymTensorCombine(alpha * strain_magnitude_filtered,
                                                  strain_filtered, 0.0, strain_filtered);
 
    // Term two: the grid-scale model formed pointwise and then test-filtered. This is
    // the filter of a product, and it is not equal to the product of the filtered
    // factors above. That inequality carries the whole information content of the
    // dynamic procedure.
    const SymTensor grid_scale = strain_product_filtered;
 
    const SymTensor model = SymTensorCombine(-2.0 * delta * delta, test_scale,
                                             2.0 * delta * delta, grid_scale);
 
    // Trace-free, so the Leonard stress's trace cannot enter the contraction through
    // whatever discrete divergence the velocity field carries.
    return SymTensorDeviator(model);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ResolveLESAveragingDirections()

PetscErrorCode ResolveLESAveragingDirections	(	UserCtx *	user,
		PetscBool	direction[3]
	)

Resolves which logical directions the dynamic coefficient is averaged over.

The resulting selector is passed to PicurvSpatialRatioAverage(), which performs the reduction; this routine only decides which directions it spans.

LES_AVERAGING_LOCAL selects none, LES_AVERAGING_GLOBAL selects all three, and LES_AVERAGING_HOMOGENEOUS uses the directions the configuration names, defaulting to the block's periodic axes when it names none. Periodicity is read from the boundary face configuration, which is the same source the periodic field synchronization uses.

Parameters

[in]	user	Block context supplying configuration and boundary faces.
[out]	direction	Per-direction averaging flags in (xi, eta, zeta) order.

Returns

PetscErrorCode 0 on success, or PETSC_ERR_ARG_WRONGSTATE when homogeneous averaging is requested but no direction can be resolved.

Resolves which logical directions the dynamic coefficient is averaged over.

Full API contract is documented with the header declaration in include/les.h.

See also

ResolveLESAveragingDirections()

Definition at line 410 of file les.c.

{
    const LESConfig *config = &user->simCtx->les_config;
    PetscInt         selected = 0;
    PetscBool        periodic[3];
 
    PetscFunctionBeginUser;
    PetscCheck(direction != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Averaging direction destination cannot be NULL.");
    LESPeriodicAxes(user->simCtx, periodic);
 
    for (PetscInt axis = 0; axis < 3; axis++) direction[axis] = PETSC_FALSE;
 
    switch (config->averaging_mode) {
    case LES_AVERAGING_LOCAL:
        // No averaging: the coefficient is whatever the single cell's contractions give.
        PetscFunctionReturn(0);
 
    case LES_AVERAGING_GLOBAL:
        for (PetscInt axis = 0; axis < 3; axis++) direction[axis] = PETSC_TRUE;
        PetscFunctionReturn(0);
 
    case LES_AVERAGING_HOMOGENEOUS:
    default:
        for (PetscInt axis = 0; axis < 3; axis++) {
            if (config->averaging_direction[axis]) {
                direction[axis] = PETSC_TRUE;
                selected++;
            }
        }
        // Naming no direction asks for the periodic axes, which is where the flow is
        // homogeneous in every case the solver can currently express.
        if (selected == 0) {
            for (PetscInt axis = 0; axis < 3; axis++) {
                if (periodic[axis]) {
                    direction[axis] = PETSC_TRUE;
                    selected++;
                }
            }
        }
        PetscCheck(selected > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                   "Homogeneous LES averaging was requested, but this block declares no "
                   "periodic axis and the configuration names no direction. Name the "
                   "homogeneous directions explicitly or select local averaging.");
        PetscFunctionReturn(0);
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

ClipModelCoefficient()

PetscReal ClipModelCoefficient	(	PetscReal	coefficient,
		const LESConfig *	config,
		PetscBool *	limited
	)

Applies the configured admissible range to one model coefficient.

The coefficient carries the sign of <L:M>; a negative value is backscatter, energy returning from the unresolved scales, not an error. LES_CLIP_NONE keeps it, leaving stability to the total-viscosity floor applied downstream.

Parameters

[in]	coefficient	Raw coefficient from the Germano contraction.
[in]	config	LES configuration supplying the mode and ceiling.
[out]	limited	Set to PETSC_TRUE when the value was modified; may be NULL.

Returns

The coefficient after limiting.

Applies the configured admissible range to one model coefficient.

Full API contract is documented with the header declaration in include/les.h.

See also

ClipModelCoefficient()

Definition at line 466 of file les.c.

{
    PetscReal result = coefficient;
 
    switch (config->clip_mode) {
    case LES_CLIP_NONE:
        // Signed coefficient retained: backscatter survives, and stability is left to
        // the total-viscosity floor rather than to the coefficient's sign.
        break;
    case LES_CLIP_CLIP_NEGATIVE:
        if (result < 0.0) result = 0.0;
        break;
    case LES_CLIP_CLAMP:
    default:
        if (result < 0.0) result = 0.0;
        if (result > config->max_cs * config->max_cs) result = config->max_cs * config->max_cs;
        break;
    }
 
    if (limited) *limited = (PetscBool)(result != coefficient);
    return result;
}

Here is the caller graph for this function:

EddyViscosityFromCoefficient()

PetscReal EddyViscosityFromCoefficient	(	PetscReal	coefficient,
		PetscReal	delta,
		PetscReal	strain_magnitude,
		PetscReal	molecular_viscosity,
		PetscReal	min_viscosity_ratio
	)

Builds the eddy viscosity from a model coefficient.

nu_t = C Delta^2 |S|, floored so that nu + nu_t >= min_viscosity_ratio * nu. The floor constrains the quantity that must stay positive for the momentum operator to remain well posed, which is why backscatter is limited here rather than by forcing the coefficient positive.

Parameters

[in]	coefficient	Model coefficient, signed.
[in]	delta	Grid filter width at the cell.
[in]	strain_magnitude	Resolved strain magnitude |S|.
[in]	molecular_viscosity	Molecular kinematic viscosity.
[in]	min_viscosity_ratio	Lower bound on total viscosity as a multiple of molecular.

Returns

The admissible eddy viscosity.

Builds the eddy viscosity from a model coefficient.

Full API contract is documented with the header declaration in include/les.h.

See also

EddyViscosityFromCoefficient()

Definition at line 497 of file les.c.

{
    PetscReal nu_t = coefficient * delta * delta * strain_magnitude;
 
    // A negative eddy viscosity is physical up to the point where it would drive the
    // total viscosity non-positive and make the momentum operator ill posed.
    const PetscReal floor = (min_viscosity_ratio - 1.0) * molecular_viscosity;
 
    if (nu_t < floor) nu_t = floor;
    return nu_t;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WALEEddyViscosity()

PetscReal WALEEddyViscosity	(	Cmpnts	dudx,
		Cmpnts	dvdx,
		Cmpnts	dwdx,
		PetscReal	delta,
		PetscReal	coefficient
	)

WALE eddy viscosity (Nicoud & Ducros 1999) from the resolved velocity gradient.

nu_t = (C_w Delta)^2 (S^d_ij S^d_ij)^(3/2) / [(S_ij S_ij)^(5/2) + (S^d_ij S^d_ij)^(5/4)], where S^d is the traceless symmetric part of the squared gradient. It vanishes in pure shear, so a laminar shear layer gets none, and it falls as the cube of wall distance approaching a wall without a damping function. It needs no coefficient field, no test filter, and no averaging, so it is valid on geometry with no homogeneous direction.

Parameters

[in]	dudx	Gradient of the x velocity component.
[in]	dvdx	Gradient of the y velocity component.
[in]	dwdx	Gradient of the z velocity component.
[in]	delta	Grid filter width at the cell.
[in]	coefficient	Model constant C_w.

Returns

The eddy viscosity, zero where the gradient carries no model invariant.

WALE eddy viscosity (Nicoud & Ducros 1999) from the resolved velocity gradient.

Full API contract is documented with the header declaration in include/les.h.

See also

WALEEddyViscosity()

Definition at line 520 of file les.c.

{
    const PetscReal g[3][3] = {{dudx.x, dudx.y, dudx.z},
                               {dvdx.x, dvdx.y, dvdx.z},
                               {dwdx.x, dwdx.y, dwdx.z}};
    PetscReal g2[3][3], sd[3][3], trace = 0.0, ss = 0.0, sdsd = 0.0, denominator;
 
    for (PetscInt a = 0; a < 3; ++a)
        for (PetscInt b = 0; b < 3; ++b) {
            g2[a][b] = 0.0;
            for (PetscInt c = 0; c < 3; ++c) g2[a][b] += g[a][c] * g[c][b];
        }
    for (PetscInt a = 0; a < 3; ++a) trace += g2[a][a];
    for (PetscInt a = 0; a < 3; ++a)
        for (PetscInt b = 0; b < 3; ++b) {
            const PetscReal s = 0.5 * (g[a][b] + g[b][a]);
            // Traceless symmetric part of the squared gradient: it vanishes in pure shear,
            // which is what lets the model switch itself off in a laminar shear layer.
            sd[a][b] = 0.5 * (g2[a][b] + g2[b][a]) - (a == b ? trace / 3.0 : 0.0);
            ss   += s * s;
            sdsd += sd[a][b] * sd[a][b];
        }
 
    denominator = PetscPowReal(ss, 2.5) + PetscPowReal(sdsd, 1.25);
    if (denominator <= LES_EPSILON * LES_EPSILON) return 0.0;
    return coefficient * coefficient * delta * delta * PetscPowReal(sdsd, 1.5) / denominator;
}

Here is the caller graph for this function:

VremanEddyViscosity()

PetscReal VremanEddyViscosity	(	Cmpnts	dudx,
		Cmpnts	dvdx,
		Cmpnts	dwdx,
		const Cmpnts	edges[3],
		PetscReal	coefficient
	)

Vreman eddy viscosity (Vreman 2004), weighting each grid direction by its own spacing.

nu_t = c sqrt(B_beta / (alpha_ij alpha_ij)), with alpha_ij = du_j/dx_i and beta_ij = sum_m d_i^(m) d_j^(m), where d^(m) is the change in velocity across the cell along grid direction m - edges[m] . grad(u). On an aligned Cartesian grid that is Vreman's Delta_m^2 alpha_mi alpha_mj; on a curvilinear grid it is the same quantity measured along the cell's own edges, so the model is resolved per direction rather than through one scalar width, and unchanged by the grid's orientation. It vanishes in pure shear and needs no coefficient field, test filter, or averaging.

Parameters

[in]	dudx	Gradient of the x velocity component.
[in]	dvdx	Gradient of the y velocity component.
[in]	dwdx	Gradient of the z velocity component.
[in]	edges	Cell edge vectors along xi, eta, zeta, from ComputeCellEdgeVectors.
[in]	coefficient	Model constant c.

Returns

The eddy viscosity, zero where B_beta vanishes.

Vreman eddy viscosity (Vreman 2004), weighting each grid direction by its own spacing.

Full API contract is documented with the header declaration in include/les.h.

See also

VremanEddyViscosity()

Definition at line 557 of file les.c.

{
    const Cmpnts gradient[3] = {dudx, dvdx, dwdx};
    PetscReal beta[3][3] = {{0.0}}, alpha_sq = 0.0, b_beta;
 
    // Vreman's beta_ij = sum_m Delta_m^2 alpha_mi alpha_mj on a Cartesian grid, where
    // Delta_m alpha_mi is the change in u_i across the cell in direction m. On a
    // curvilinear grid that change is edge_m . grad(u_i), so the same sum over the
    // cell's own directions reproduces the Cartesian form on an aligned grid and does
    // not depend on the grid's orientation elsewhere.
    for (PetscInt m = 0; m < 3; ++m) {
        PetscReal change[3];
        for (PetscInt a = 0; a < 3; ++a)
            change[a] = edges[m].x * gradient[a].x + edges[m].y * gradient[a].y + edges[m].z * gradient[a].z;
        for (PetscInt a = 0; a < 3; ++a)
            for (PetscInt b = 0; b < 3; ++b) beta[a][b] += change[a] * change[b];
    }
    for (PetscInt a = 0; a < 3; ++a)
        alpha_sq += gradient[a].x * gradient[a].x + gradient[a].y * gradient[a].y + gradient[a].z * gradient[a].z;
 
    b_beta = beta[0][0] * beta[1][1] - beta[0][1] * beta[0][1]
           + beta[0][0] * beta[2][2] - beta[0][2] * beta[0][2]
           + beta[1][1] * beta[2][2] - beta[1][2] * beta[1][2];
 
    // B_beta is non-negative in exact arithmetic; rounding can leave it a hair below.
    if (alpha_sq <= LES_EPSILON * LES_EPSILON || b_beta <= 0.0) return 0.0;
    return coefficient * PetscSqrtReal(b_beta / alpha_sq);
}

Here is the caller graph for this function:

SubgridKineticEnergy()

PetscReal SubgridKineticEnergy	(	PetscReal	yoshizawa_ci,
		PetscReal	delta,
		PetscReal	strain_magnitude
	)

Returns the modelled subgrid kinetic energy at a cell.

Uses the Yoshizawa relation k_sgs = 2 C_I Delta^2 |S|^2. Reported as a diagnostic: compared against the resolved turbulent kinetic energy from the field-statistics pipeline it measures how much of the energy the grid resolves.

Parameters

[in]	yoshizawa_ci	Yoshizawa constant C_I.
[in]	delta	Grid filter width at the cell.
[in]	strain_magnitude	Resolved strain magnitude |S|.

Returns

The modelled subgrid kinetic energy.

Returns the modelled subgrid kinetic energy at a cell.

Full API contract is documented with the header declaration in include/les.h.

See also

SubgridKineticEnergy()

Definition at line 595 of file les.c.

{
    return 2.0 * yoshizawa_ci * delta * delta * strain_magnitude * strain_magnitude;
}

Here is the caller graph for this function:

ComputeSmagorinskyConstant()

PetscErrorCode ComputeSmagorinskyConstant

(

UserCtx *

user

)

Computes the dynamic Smagorinsky coefficient field for one block.

Runs the Germano-Lilly procedure: test-filters the resolved field, forms the Leonard stress and the model tensor, contracts them, averages the two contractions over the configured directions, limits the quotient, and stores it in UserCtx::CS.

The stored quantity is the coefficient C that multiplies Delta^2 |S| in the eddy viscosity, that is, Cs^2 rather than Cs. It is signed under LES_CLIP_NONE. Reported values are converted to Cs at the diagnostic boundary, where the familiar target near 0.16 applies.

Only DYNAMIC_SMAGORINSKY reaches this routine; the constant model prescribes its coefficient from configuration and allocates no field. Callers honour LESConfig::dynamic_frequency, and the caller must have refreshed UserCtx::lUcat and its periodic images beforehand.

Parameters

[in,out]

user

The user context for a single computational block.

Returns

PetscErrorCode 0 on success.

Computes the dynamic Smagorinsky coefficient field for one block.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/les.h.

See also

ComputeSmagorinskyConstant()

Definition at line 715 of file les.c.

{
    SimCtx          *simCtx = user->simCtx;
    const LESConfig *config = &simCtx->les_config;
    DM                da = user->da, fda = user->fda;
    SpatialTargetPlan plan;
    PetscBool         average_direction[3];
    PetscReal        alpha;
 
    Vec lStrainDiag, lStrainOff, lStrainAbs, lLM, lMM, lCsq;
    Vec lWallMask = NULL;
 
    Cmpnts    ***ucat = NULL, ***strain_diag = NULL, ***strain_off = NULL;
    Cmpnts    ***csi = NULL, ***eta = NULL, ***zet = NULL;
    PetscReal ***strain_abs = NULL, ***nvert = NULL, ***aj = NULL;
    PetscReal ***lm = NULL, ***mm = NULL, ***csq = NULL, ***cs_out = NULL;
 
    LESDiagnosticsState diagnostics = {0.0, 0.0, 0.0, 0.0, 0.0, PETSC_TRUE};
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    PetscCheck(simCtx->les == DYNAMIC_SMAGORINSKY, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "The dynamic procedure runs only for the dynamic Smagorinsky model. The "
               "constant model prescribes its coefficient from configuration and holds "
               "no coefficient field.");
 
    // The dynamic procedure needs a developed field to sample. Starting from rest it
    // has nothing to measure, so the first steps run with no subgrid dissipation.
    if (simCtx->step < 2 && simCtx->StartStep == 0) {
        PetscCall(VecSet(user->CS, 0.0));
        PetscCall(FinalizeSmagorinskyConstantField(user));
        user->les_diagnostics = (LESDiagnosticsState){0.0, 0.0, 0.0, 0.0, 0.0, PETSC_FALSE};
        LOG_ALLOW(GLOBAL, LOG_DEBUG,
                  "Holding the dynamic coefficient at zero for the opening steps (step=%d).\n",
                  simCtx->step);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    PetscCall(LESResolveDomain(user, &plan));
    PetscCall(ResolveLESAveragingDirections(user, average_direction));
 
    // alpha is the squared ratio of the two filter widths and sets the relative weight
    // of the test-scale and grid-scale terms in the model tensor. The widths are the
    // cube roots of the three directional widths, and the configured ratio applies to
    // each direction the kernel filters. The box filters all three, so alpha is the
    // ratio squared. Simpson filters only i and k and leaves j at the grid width, so
    // its widths differ by ratio^(2/3) and alpha is ratio^(4/3); squaring the ratio
    // there halved the dynamic coefficient on HOM02 (0.091 against the box's 0.194).
    if (config->test_filter_kernel == LES_TEST_FILTER_SIMPSON_IK) {
        alpha = PetscPowReal(config->test_filter_width_ratio, 4.0 / 3.0);
    } else {
        alpha = config->test_filter_width_ratio * config->test_filter_width_ratio;
    }
 
    PetscCall(VecDuplicate(user->lUcat, &lStrainDiag));
    PetscCall(VecDuplicate(user->lUcat, &lStrainOff));
    PetscCall(VecDuplicate(user->lNvert, &lStrainAbs));
    PetscCall(VecDuplicate(user->lNvert, &lLM));
    PetscCall(VecDuplicate(user->lNvert, &lMM));
    PetscCall(VecDuplicate(user->lNvert, &lCsq));
    PetscCall(VecSet(lStrainDiag, 0.0));
    PetscCall(VecSet(lStrainOff, 0.0));
    PetscCall(VecSet(lStrainAbs, 0.0));
    PetscCall(VecSet(lLM, 0.0));
    PetscCall(VecSet(lMM, 0.0));
    PetscCall(VecSet(lCsq, 0.0));
    /* A wall model overwrites the first interior cell, so those cells carry an imposed
       profile rather than a measurement and are kept out of the contraction sums. */
    if (simCtx->wallfunction) {
        PetscCall(VecDuplicate(user->lNvert, &lWallMask));
        PetscCall(BuildWallModelExclusionMask(user, lWallMask));
    }
 
    PetscCall(DMDAVecGetArray(fda, user->lUcat, &ucat));
 
    // --- 1. Strain rate everywhere the filter stencil will look -----------------
    PetscCall(ComputeStrainRateField(user, ucat, lStrainDiag, lStrainOff, lStrainAbs));
 
    PetscCall(DMDAVecGetArrayRead(fda, lStrainDiag, &strain_diag));
    PetscCall(DMDAVecGetArrayRead(fda, lStrainOff, &strain_off));
    PetscCall(DMDAVecGetArrayRead(da, lStrainAbs, &strain_abs));
    PetscCall(DMDAVecGetArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecGetArrayRead(da, user->lAj, &aj));
    PetscCall(DMDAVecGetArray(da, lLM, &lm));
    PetscCall(DMDAVecGetArray(da, lMM, &mm));
 
    // --- 2. Germano contractions at every fluid cell ----------------------------
    for (PetscInt k = plan.start[2]; k < plan.end[2]; k++)
    for (PetscInt j = plan.start[1]; j < plan.end[1]; j++)
    for (PetscInt i = plan.start[0]; i < plan.end[0]; i++) {
        double    velocity_x[3][3][3], velocity_y[3][3][3], velocity_z[3][3][3];
        double    strain_abs_stencil[3][3][3], weights[3][3][3];
        SymTensor strain_stencil[3][3][3], product_stencil[3][3][3], velocity_product[3][3][3];
        Cmpnts    velocity_filtered;
        SymTensor strain_filtered, product_filtered, velocity_product_filtered;
        SymTensor leonard, model;
        PetscReal strain_abs_filtered, delta, volume;
 
        if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) {
            lm[k][j][i] = 0.0;
            mm[k][j][i] = 0.0;
            continue;
        }
 
        // --- 2a. Gather the 3x3x3 stencil -------------------------------------
        for (PetscInt r = -1; r <= 1; r++)
        for (PetscInt q = -1; q <= 1; q++)
        for (PetscInt p = -1; p <= 1; p++) {
            const PetscInt  R = r + 1, Q = q + 1, P = p + 1;
            const PetscInt  KK = k + r, JJ = j + q, II = i + p;
            const SymTensor strain = LESStrainAt(strain_diag[KK][JJ][II], strain_off[KK][JJ][II]);
            const PetscReal strain_abs_here = strain_abs[KK][JJ][II];
 
            velocity_x[R][Q][P] = ucat[KK][JJ][II].x;
            velocity_y[R][Q][P] = ucat[KK][JJ][II].y;
            velocity_z[R][Q][P] = ucat[KK][JJ][II].z;
 
            strain_stencil[R][Q][P]     = strain;
            strain_abs_stencil[R][Q][P] = strain_abs_here;
 
            // The grid-scale model formed pointwise, before any test filtering. Test
            // filtering this product is the term the model tensor cannot do without,
            // and it is not the same as multiplying the two filtered factors.
            product_stencil[R][Q][P] = SymTensorCombine(strain_abs_here, strain, 0.0, strain);
 
            velocity_product[R][Q][P] = SymTensorSelfOuter(ucat[KK][JJ][II]);
            weights[R][Q][P]          = aj[KK][JJ][II];
        }
 
        // --- 2b. Test-filter every quantity the identity needs ----------------
        velocity_filtered.x = ApplyLESTestFilter(config->test_filter_kernel, velocity_x, weights);
        velocity_filtered.y = ApplyLESTestFilter(config->test_filter_kernel, velocity_y, weights);
        velocity_filtered.z = ApplyLESTestFilter(config->test_filter_kernel, velocity_z, weights);
        strain_abs_filtered = ApplyLESTestFilter(config->test_filter_kernel, strain_abs_stencil, weights);
        PetscCall(ApplyLESTestFilterSymTensor(config->test_filter_kernel, strain_stencil, weights, &strain_filtered));
        PetscCall(ApplyLESTestFilterSymTensor(config->test_filter_kernel, product_stencil, weights, &product_filtered));
        PetscCall(ApplyLESTestFilterSymTensor(config->test_filter_kernel, velocity_product, weights, &velocity_product_filtered));
 
        // --- 2c. Leonard stress and model tensor -------------------------------
        leonard = LeonardStress(velocity_filtered, velocity_product_filtered);
 
        PetscCall(ComputeCellFilterWidth(config->filter_width_model, aj[k][j][i],
                                         csi[k][j][i], eta[k][j][i], zet[k][j][i], &delta));
 
        model = GermanoModelTensor(delta, alpha, strain_abs_filtered, strain_filtered, product_filtered);
 
        // --- 2d. Lilly's least-squares contraction -----------------------------
        // Stored pre-multiplied by cell volume. The average that follows is over
        // volume, not over cell count, so that a refined region does not dominate it
        // merely by holding more cells; carrying the weight here keeps the shared
        // averaging kernel free of any weighting policy of its own.
        volume = 1.0 / aj[k][j][i];
        lm[k][j][i] = volume * SymTensorContract(leonard, model);
        mm[k][j][i] = volume * SymTensorNormSq(model);
    }
 
    PetscCall(DMDAVecRestoreArrayRead(fda, lStrainDiag, &strain_diag));
    PetscCall(DMDAVecRestoreArrayRead(fda, lStrainOff, &strain_off));
    PetscCall(DMDAVecRestoreArray(da, lLM, &lm));
    PetscCall(DMDAVecRestoreArray(da, lMM, &mm));
    PetscCall(DMDAVecRestoreArray(fda, user->lUcat, &ucat));
 
    // --- 3. Average the two contractions, then divide ---------------------------
    // Both sums are formed before the division. Averaging the pointwise quotients
    // instead would not be the least-squares coefficient Lilly's closure defines.
    PetscCall(PicurvSpatialRatioAverage(user, &plan, lLM, lMM, lWallMask, average_direction,
                                        PetscObjectComm((PetscObject)da), lCsq, NULL));
 
    // --- 4. Limit the coefficient and record what limiting cost -----------------
    PetscCall(DMDAVecGetArrayRead(da, lLM, &lm));
    PetscCall(DMDAVecGetArrayRead(da, lMM, &mm));
    PetscCall(DMDAVecGetArrayRead(da, lCsq, &csq));
    PetscCall(DMDAVecGetArray(da, user->CS, &cs_out));
 
    for (PetscInt k = plan.start[2]; k < plan.end[2]; k++)
    for (PetscInt j = plan.start[1]; j < plan.end[1]; j++)
    for (PetscInt i = plan.start[0]; i < plan.end[0]; i++) {
        PetscBool limited = PETSC_FALSE;
        PetscReal weight, raw;
 
        if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) {
            cs_out[k][j][i] = 0.0;
            continue;
        }
 
        weight = 1.0 / aj[k][j][i];
        raw    = csq[k][j][i];
 
        cs_out[k][j][i] = ClipModelCoefficient(raw, config, &limited);
 
        // The two contraction fields already carry the cell volume from the loop
        // above, so they are summed as they stand rather than weighted twice.
        diagnostics.contraction_lm     += lm[k][j][i];
        diagnostics.contraction_mm     += mm[k][j][i];
        diagnostics.fluid_volume       += weight;
        if (raw < 0.0) diagnostics.backscatter_volume += weight;
        if (limited)   diagnostics.limited_volume     += weight;
    }
 
    PetscCall(DMDAVecRestoreArrayRead(da, lLM, &lm));
    PetscCall(DMDAVecRestoreArrayRead(da, lMM, &mm));
    PetscCall(DMDAVecRestoreArrayRead(da, lCsq, &csq));
    PetscCall(DMDAVecRestoreArray(da, user->CS, &cs_out));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lAj, &aj));
 
    PetscCall(VecDestroy(&lStrainDiag));
    PetscCall(VecDestroy(&lStrainOff));
    PetscCall(VecDestroy(&lStrainAbs));
    PetscCall(VecDestroy(&lLM));
    PetscCall(VecDestroy(&lMM));
    PetscCall(VecDestroy(&lCsq));
    if (lWallMask) PetscCall(VecDestroy(&lWallMask));
 
    user->les_diagnostics = diagnostics;
    PetscCall(FinalizeSmagorinskyConstantField(user));
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeEddyViscosityLES()

PetscErrorCode ComputeEddyViscosityLES

(

UserCtx *

user

)

Computes the turbulent eddy viscosity for one block.

Evaluates nu_t = C Delta^2 |S| at every fluid cell and stores it in UserCtx::Nu_t. The coefficient comes from UserCtx::lCs for the dynamic model and directly from LESConfig::constant_cs for the constant model, which therefore needs no coefficient field at all.

Parameters

[in,out]

user

The user context for a single computational block.

Returns

PetscErrorCode 0 on success.

Computes the turbulent eddy viscosity for one block.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/les.h.

See also

ComputeEddyViscosityLES()

Definition at line 952 of file les.c.

{
    SimCtx          *simCtx = user->simCtx;
    const LESConfig *config = &simCtx->les_config;
    DM                da = user->da, fda = user->fda;
    SpatialTargetPlan plan;
    const PetscReal  molecular_viscosity = 1.0 / simCtx->ren;
 
    // The constant model's coefficient is a number from the configuration, so it needs
    // no field, no ghost exchange, and no periodic synchronization.
    const PetscBool  coefficient_is_field = (PetscBool)(simCtx->les == DYNAMIC_SMAGORINSKY);
    const PetscReal  prescribed_coefficient = config->constant_cs * config->constant_cs;
 
    Cmpnts    ***ucat = NULL, ***csi = NULL, ***eta = NULL, ***zet = NULL;
    PetscReal ***coefficient = NULL, ***nu_t = NULL, ***nvert = NULL, ***aj = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    PetscCall(LESResolveDomain(user, &plan));
 
    PetscCall(DMDAVecGetArrayRead(fda, user->lUcat, &ucat));
    PetscCall(DMDAVecGetArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArray(da, user->Nu_t, &nu_t));
    PetscCall(DMDAVecGetArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecGetArrayRead(da, user->lAj, &aj));
    if (coefficient_is_field) PetscCall(DMDAVecGetArrayRead(da, user->lCs, &coefficient));
 
    for (PetscInt k = plan.start[2]; k < plan.end[2]; k++)
    for (PetscInt j = plan.start[1]; j < plan.end[1]; j++)
    for (PetscInt i = plan.start[0]; i < plan.end[0]; i++) {
        Cmpnts    dudx, dvdx, dwdx;
        PetscReal strain_magnitude, delta, model_coefficient;
 
        if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) {
            nu_t[k][j][i] = 0.0;
            continue;
        }
 
        PetscCall(ComputeVectorFieldDerivatives(user, i, j, k, (Cmpnts ***)ucat, &dudx, &dvdx, &dwdx));
 
        if (simCtx->les == VREMAN) {
            // Resolved per direction from the cell's own edges; no scalar width enters.
            Cmpnts edges[3];
            PetscCall(ComputeCellEdgeVectors(aj[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i], edges));
            nu_t[k][j][i] = VremanEddyViscosity(dudx, dvdx, dwdx, edges, config->vreman_coefficient);
            continue;
        }
 
        PetscCall(ComputeCellFilterWidth(config->filter_width_model, aj[k][j][i],
                                         csi[k][j][i], eta[k][j][i], zet[k][j][i], &delta));
        if (simCtx->les == WALE) {
            nu_t[k][j][i] = WALEEddyViscosity(dudx, dvdx, dwdx, delta, config->wale_coefficient);
            continue;
        }
 
        PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, NULL, &strain_magnitude));
        model_coefficient = coefficient_is_field ? coefficient[k][j][i] : prescribed_coefficient;
 
        nu_t[k][j][i] = EddyViscosityFromCoefficient(model_coefficient, delta, strain_magnitude,
                                                     molecular_viscosity, config->min_viscosity_ratio);
    }
 
    if (coefficient_is_field) PetscCall(DMDAVecRestoreArrayRead(da, user->lCs, &coefficient));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lUcat, &ucat));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArray(da, user->Nu_t, &nu_t));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lAj, &aj));
 
    {
        const FieldId periodic_fields[] = {FIELD_ID_NU_T};
 
        PetscCall(SynchronizePeriodicCellFields(user, 1, periodic_fields));
        PetscCall(UpdateLocalGhosts(user, FIELD_ID_NU_T));
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LogLESDiagnostics()

PetscErrorCode LogLESDiagnostics

(

UserCtx *

user

)

Appends one row of LES coefficient statistics to the run's log directory.

Writes les_coefficient.csv: the effective coefficient reported as Cs, its spatial spread, eddy-viscosity levels, the modelled subgrid energy, and the fractions of the domain that were backscattering or limited before this step's clipping. Those last two describe the pre-clipping state, which no stored field preserves.

Values are instantaneous volume statistics, not time averages; window averages of the same model's fields come from the field-statistics pipeline instead.

Parameters

[in]

user

The user context for a single computational block.

Returns

PetscErrorCode 0 on success.

Appends one row of LES coefficient statistics to the run's log directory.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/les.h.

See also

LogLESDiagnostics()

Definition at line 1045 of file les.c.

{
    SimCtx          *simCtx = user->simCtx;
    const LESConfig *config = &simCtx->les_config;
    DM                da = user->da, fda = user->fda;
    SpatialTargetPlan plan;
    MPI_Comm         comm;
 
    const PetscBool dynamic = (PetscBool)(simCtx->les == DYNAMIC_SMAGORINSKY);
    // Vreman and WALE carry no Smagorinsky coefficient, so the coefficient columns are
    // written as nan rather than as a number that would read like a measured Cs.
    const PetscBool has_coefficient =
        (PetscBool)(simCtx->les == DYNAMIC_SMAGORINSKY || simCtx->les == CONSTANT_SMAGORINSKY);
 
    // Index 0..4: fluid volume, coefficient first and second moments, nu_t volume sum,
    // subgrid-energy volume sum. Index 5..6: the pre-clip fractions. Reduced together
    // so the diagnostic costs one collective rather than seven.
    PetscReal local_sum[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    PetscReal global_sum[7];
    PetscReal local_max[2] = {0.0, 0.0}, global_max[2];
    PetscReal local_min_coefficient = PETSC_MAX_REAL, global_min_coefficient;
    PetscReal contraction[2], global_contraction[2];
 
    Cmpnts    ***ucat = NULL, ***csi = NULL, ***eta = NULL, ***zet = NULL;
    PetscReal ***coefficient = NULL, ***nu_t = NULL, ***nvert = NULL, ***aj = NULL;
 
    PetscFunctionBeginUser;
 
    if (!config->diagnostics_enabled) PetscFunctionReturn(0);
    if (config->diagnostics_cadence <= 0) PetscFunctionReturn(0);
    if (simCtx->step % config->diagnostics_cadence != 0) PetscFunctionReturn(0);
 
    PetscCall(LESResolveDomain(user, &plan));
    PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
 
    PetscCall(DMDAVecGetArrayRead(fda, user->lUcat, &ucat));
    PetscCall(DMDAVecGetArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArrayRead(da, user->lNu_t, &nu_t));
    PetscCall(DMDAVecGetArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecGetArrayRead(da, user->lAj, &aj));
    if (dynamic) PetscCall(DMDAVecGetArrayRead(da, user->lCs, &coefficient));
 
    for (PetscInt k = plan.start[2]; k < plan.end[2]; k++)
    for (PetscInt j = plan.start[1]; j < plan.end[1]; j++)
    for (PetscInt i = plan.start[0]; i < plan.end[0]; i++) {
        Cmpnts    dudx, dvdx, dwdx;
        PetscReal strain_magnitude, delta, weight, value;
 
        if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) continue;
 
        weight = 1.0 / aj[k][j][i];
        value  = dynamic ? coefficient[k][j][i]
                         : config->constant_cs * config->constant_cs;
 
        PetscCall(ComputeVectorFieldDerivatives(user, i, j, k, (Cmpnts ***)ucat, &dudx, &dvdx, &dwdx));
        PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, NULL, &strain_magnitude));
        PetscCall(ComputeCellFilterWidth(config->filter_width_model, aj[k][j][i],
                                         csi[k][j][i], eta[k][j][i], zet[k][j][i], &delta));
 
        local_sum[0] += weight;
        local_sum[3] += weight * nu_t[k][j][i];
        local_sum[4] += weight * SubgridKineticEnergy(config->yoshizawa_ci, delta, strain_magnitude);
        local_max[1] = PetscMax(local_max[1], nu_t[k][j][i]);
        if (has_coefficient) {
            local_sum[1] += weight * value;
            local_sum[2] += weight * value * value;
            local_max[0] = PetscMax(local_max[0], value);
            local_min_coefficient = PetscMin(local_min_coefficient, value);
        }
    }
 
    if (dynamic) PetscCall(DMDAVecRestoreArrayRead(da, user->lCs, &coefficient));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lUcat, &ucat));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lCsi, &csi));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lNu_t, &nu_t));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lNvert, &nvert));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lAj, &aj));
 
    local_sum[5] = user->les_diagnostics.valid ? user->les_diagnostics.backscatter_volume : 0.0;
    local_sum[6] = user->les_diagnostics.valid ? user->les_diagnostics.limited_volume : 0.0;
    contraction[0] = user->les_diagnostics.valid ? user->les_diagnostics.contraction_lm : 0.0;
    contraction[1] = user->les_diagnostics.valid ? user->les_diagnostics.contraction_mm : 0.0;
 
    PetscCallMPI(MPI_Allreduce(local_sum, global_sum, 7, MPIU_REAL, MPI_SUM, comm));
    PetscCallMPI(MPI_Allreduce(local_max, global_max, 2, MPIU_REAL, MPI_MAX, comm));
    PetscCallMPI(MPI_Allreduce(&local_min_coefficient, &global_min_coefficient, 1, MPIU_REAL, MPI_MIN, comm));
    PetscCallMPI(MPI_Allreduce(contraction, global_contraction, 2, MPIU_REAL, MPI_SUM, comm));
 
    if (simCtx->rank == 0) {
        const PetscReal volume        = (global_sum[0] > 0.0) ? global_sum[0] : 1.0;
        const PetscReal mean          = global_sum[1] / volume;
        const PetscReal mean_square   = global_sum[2] / volume;
        const PetscReal variance      = PetscMax(mean_square - mean * mean, 0.0);
        const PetscReal molecular     = 1.0 / simCtx->ren;
        const PetscReal effective     = (global_contraction[1] > LES_EPSILON)
                                            ? (global_contraction[0] / global_contraction[1]) : 0.0;
        FILE *file = NULL;
 
        // Reported as Cs rather than as the stored coefficient, because that is the
        // number with a recognised value: near 0.16 for decaying isotropic turbulence.
        // A negative average means the box is backscattering on balance, and is written
        // as a negative Cs rather than hidden behind a square root.
        const PetscReal cs_effective = (effective >= 0.0) ? PetscSqrtReal(effective)
                                                          : -PetscSqrtReal(-effective);
        const PetscReal cs_mean      = (mean >= 0.0) ? PetscSqrtReal(mean) : -PetscSqrtReal(-mean);
 
        PetscCall(PicurvOpenDiagnosticsCsv(simCtx, "les_coefficient.csv",
                                           "step,time,cs_effective,cs_mean,coefficient_mean,"
                                           "coefficient_rms,coefficient_min,coefficient_max,"
                                           "nu_t_mean,nu_t_max,nu_t_over_nu_mean,k_sgs_mean,"
                                           "backscatter_fraction,limited_fraction,physical_time", &file));
        PetscReal physical_time = 0.0;
        PetscCall(PicurvPhysicalTime(simCtx, simCtx->ti, &physical_time));
        fprintf(file,
                "%d,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e\n",
                (int)simCtx->step, (double)simCtx->ti,
                has_coefficient ? (double)cs_effective : (double)NAN,
                has_coefficient ? (double)cs_mean : (double)NAN,
                has_coefficient ? (double)mean : (double)NAN,
                has_coefficient ? (double)PetscSqrtReal(variance) : (double)NAN,
                has_coefficient ? (double)global_min_coefficient : (double)NAN,
                has_coefficient ? (double)global_max[0] : (double)NAN,
                (double)(global_sum[3] / volume), (double)global_max[1],
                (double)(global_sum[3] / volume / molecular),
                (double)(global_sum[4] / volume),
                (double)(global_sum[5] / volume), (double)(global_sum[6] / volume),
                (double)physical_time);
        PetscCheck(fclose(file) == 0, PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE,
                   "Unable to close the LES diagnostics file.");
 
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "  LES coefficient: Cs_effective=%.4f, nu_t/nu (mean)=%.3e, backscatter=%.1f%%\n",
                  (double)cs_effective, (double)(global_sum[3] / volume / molecular),
                  (double)(100.0 * global_sum[5] / volume));
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: