les.c File Reference

Implements the Large Eddy Simulation (LES) subgrid-scale closure. More...

#include "les.h"
#include "statistics_target.h"

Include dependency graph for les.c:

Go to the source code of this file.

Macros
#define	__FUNCT__   "SymTensorSelfOuter"
#define	__FUNCT__   "SymTensorCombine"
#define	__FUNCT__   "SymTensorTrace"
#define	__FUNCT__   "SymTensorDeviator"
#define	__FUNCT__   "SymTensorContract"
#define	__FUNCT__   "SymTensorTimesVector"
#define	__FUNCT__   "SymTensorNormSq"
#define	__FUNCT__   "StrainRateFromGradients"
#define	__FUNCT__   "ComputeCellFilterWidth"
#define	__FUNCT__   "LeonardStress"
#define	__FUNCT__   "GermanoModelTensor"
#define	__FUNCT__   "ResolveLESAveragingDirections"
#define	__FUNCT__   "ClipModelCoefficient"
#define	__FUNCT__   "EddyViscosityFromCoefficient"
#define	__FUNCT__   "WALEEddyViscosity"
#define	__FUNCT__   "VremanEddyViscosity"
#define	__FUNCT__   "SubgridKineticEnergy"
#define	__FUNCT__   "ComputeSmagorinskyConstant"
#define	__FUNCT__   "ComputeEddyViscosityLES"
#define	__FUNCT__   "LogLESDiagnostics"

Functions
static PetscErrorCode	FinalizeSmagorinskyConstantField (UserCtx *user)
	Synchronizes the completed Smagorinsky coefficient field.
static PetscErrorCode	LESResolveDomain (UserCtx *user, SpatialTargetPlan *plan)
	Resolves the iteration domain the closure computes over.
static PetscErrorCode	BuildWallModelExclusionMask (UserCtx *user, Vec inclusion)
	Marks the near-wall cells a wall model overwrites, so the procedure skips them.
static void	LESPeriodicAxes (const SimCtx *simCtx, PetscBool periodic[3])
	Reports the block's per-axis periodicity in (xi, eta, zeta) order.
SymTensor	SymTensorSelfOuter (Cmpnts v)
	Implementation of SymTensorSelfOuter().
SymTensor	SymTensorCombine (PetscReal a, SymTensor x, PetscReal b, SymTensor y)
	Implementation of SymTensorCombine().
PetscReal	SymTensorTrace (SymTensor t)
	Implementation of SymTensorTrace().
SymTensor	SymTensorDeviator (SymTensor t)
	Implementation of SymTensorDeviator().
PetscReal	SymTensorContract (SymTensor a, SymTensor b)
	Implementation of SymTensorContract().
Cmpnts	SymTensorTimesVector (SymTensor t, Cmpnts v)
	Implementation of SymTensorTimesVector().
PetscReal	SymTensorNormSq (SymTensor t)
	Implementation of SymTensorNormSq().
PetscErrorCode	StrainRateFromGradients (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, SymTensor *strain, PetscReal *magnitude)
	Implementation of StrainRateFromGradients().
PetscErrorCode	ComputeCellFilterWidth (LESFilterWidthModel model, PetscReal aj, Cmpnts csi, Cmpnts eta, Cmpnts zet, PetscReal *delta)
	Implementation of ComputeCellFilterWidth().
SymTensor	LeonardStress (Cmpnts velocity_filtered, SymTensor velocity_product_filtered)
	Implementation of LeonardStress().
SymTensor	GermanoModelTensor (PetscReal delta, PetscReal alpha, PetscReal strain_magnitude_filtered, SymTensor strain_filtered, SymTensor strain_product_filtered)
	Implementation of GermanoModelTensor().
PetscErrorCode	ResolveLESAveragingDirections (UserCtx *user, PetscBool direction[3])
	Implementation of ResolveLESAveragingDirections().
PetscReal	ClipModelCoefficient (PetscReal coefficient, const LESConfig *config, PetscBool *limited)
	Implementation of ClipModelCoefficient().
PetscReal	EddyViscosityFromCoefficient (PetscReal coefficient, PetscReal delta, PetscReal strain_magnitude, PetscReal molecular_viscosity, PetscReal min_viscosity_ratio)
	Implementation of EddyViscosityFromCoefficient().
PetscReal	WALEEddyViscosity (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, PetscReal delta, PetscReal coefficient)
	Implementation of WALEEddyViscosity().
PetscReal	VremanEddyViscosity (Cmpnts dudx, Cmpnts dvdx, Cmpnts dwdx, const Cmpnts edges[3], PetscReal coefficient)
	Implementation of VremanEddyViscosity().
PetscReal	SubgridKineticEnergy (PetscReal yoshizawa_ci, PetscReal delta, PetscReal strain_magnitude)
	Implementation of SubgridKineticEnergy().
static PetscInt	LESStrainSampleIndex (PetscInt index, PetscInt extent, PetscBool periodic)
	Maps a cell index to the index whose strain rate represents it.
static PetscErrorCode	ComputeStrainRateField (UserCtx *user, Cmpnts ***ucat, Vec strain_diagonal, Vec strain_offdiagonal, Vec strain_magnitude)
	Fills the strain-rate scratch fields across the computed range and its halo.
static SymTensor	LESStrainAt (Cmpnts diagonal, Cmpnts offdiagonal)
	Rebuilds a symmetric strain tensor from its two scratch storage vectors.
PetscErrorCode	ComputeSmagorinskyConstant (UserCtx *user)
	Implementation of ComputeSmagorinskyConstant().
PetscErrorCode	ComputeEddyViscosityLES (UserCtx *user)
	Implementation of ComputeEddyViscosityLES().
PetscErrorCode	LogLESDiagnostics (UserCtx *user)
	Implementation of LogLESDiagnostics().

Variables
static const double	LES_EPSILON = 1.0e-12

Detailed Description

Implements the Large Eddy Simulation (LES) subgrid-scale closure.

The closure is built from small kernels that each do one thing, so the physics is readable at the call sites and testable without a running solve:

• symmetric-tensor algebra (SymTensor*), which knows nothing about turbulence;
• strain rate and filter width, shared by every eddy-viscosity model;
• the Germano identity pieces, LeonardStress() and GermanoModelTensor();
• the coefficient averaging, which is PicurvSpatialRatioAverage() from the shared spatial-target module: it turns the local dynamic model into the plane- or volume-averaged one without a separate code path per mode;
• limiting and eddy-viscosity assembly.

Two routines drive them. ComputeSmagorinskyConstant() runs the Germano-Lilly procedure for the dynamic model, and ComputeEddyViscosityLES() turns whichever coefficient is in force into the eddy viscosity the momentum equations add to the molecular one. LogLESDiagnostics() reports the coefficient's time history.

The coefficient stored in UserCtx::CS is C, the factor multiplying Delta^2 |S|, which is Cs^2 in the classical notation and is signed when backscatter is admitted. Conversion to the familiar Cs happens only where a human reads the number.

Definition in file les.c.

Macro Definition Documentation

__FUNCT__ [1/20]

#define __FUNCT__   "SymTensorSelfOuter"

Definition at line 137 of file les.c.

__FUNCT__ [2/20]

#define __FUNCT__   "SymTensorCombine"

Definition at line 137 of file les.c.

__FUNCT__ [3/20]

#define __FUNCT__   "SymTensorTrace"

Definition at line 137 of file les.c.

__FUNCT__ [4/20]

#define __FUNCT__   "SymTensorDeviator"

Definition at line 137 of file les.c.

__FUNCT__ [5/20]

#define __FUNCT__   "SymTensorContract"

Definition at line 137 of file les.c.

__FUNCT__ [6/20]

#define __FUNCT__   "SymTensorTimesVector"

Definition at line 137 of file les.c.

__FUNCT__ [7/20]

#define __FUNCT__   "SymTensorNormSq"

Definition at line 137 of file les.c.

__FUNCT__ [8/20]

#define __FUNCT__   "StrainRateFromGradients"

Definition at line 137 of file les.c.

__FUNCT__ [9/20]

#define __FUNCT__   "ComputeCellFilterWidth"

Definition at line 137 of file les.c.

__FUNCT__ [10/20]

#define __FUNCT__   "LeonardStress"

Definition at line 137 of file les.c.

__FUNCT__ [11/20]

#define __FUNCT__   "GermanoModelTensor"

Definition at line 137 of file les.c.

__FUNCT__ [12/20]

#define __FUNCT__   "ResolveLESAveragingDirections"

Definition at line 137 of file les.c.

__FUNCT__ [13/20]

#define __FUNCT__   "ClipModelCoefficient"

Definition at line 137 of file les.c.

__FUNCT__ [14/20]

#define __FUNCT__   "EddyViscosityFromCoefficient"

Definition at line 137 of file les.c.

__FUNCT__ [15/20]

#define __FUNCT__   "WALEEddyViscosity"

Definition at line 137 of file les.c.

__FUNCT__ [16/20]

#define __FUNCT__   "VremanEddyViscosity"

Definition at line 137 of file les.c.

__FUNCT__ [17/20]

#define __FUNCT__   "SubgridKineticEnergy"

Definition at line 137 of file les.c.

__FUNCT__ [18/20]

#define __FUNCT__   "ComputeSmagorinskyConstant"

Definition at line 137 of file les.c.

__FUNCT__ [19/20]

#define __FUNCT__   "ComputeEddyViscosityLES"

Definition at line 137 of file les.c.

__FUNCT__ [20/20]

#define __FUNCT__   "LogLESDiagnostics"

Definition at line 137 of file les.c.

Function Documentation

FinalizeSmagorinskyConstantField()

static PetscErrorCode FinalizeSmagorinskyConstantField ( UserCtx *  user )

static

Synchronizes the completed Smagorinsky coefficient field.

Definition at line 40 of file les.c.

{
    const FieldId fields[] = {FIELD_ID_CS};
 
    PetscFunctionBeginUser;
    PetscCall(SynchronizePeriodicCellFields(user, 1, fields));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_CS));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LESResolveDomain()

static PetscErrorCode LESResolveDomain ( UserCtx *  user, SpatialTargetPlan *  plan  )

static

Resolves the iteration domain the closure computes over.

Delegates to the shared spatial target so the closure, the statistics pipeline, and anything else looping over cell-centred data agree on which indices are physical. The domain excludes PETSc halo storage and the solver's boundary, dummy, and duplicate-periodic indices; the last of those matters here, because counting a duplicate plane in a spatial average would weight a real cell twice.

The plan is resolved for CS as a representative cell-centred field. The bounds depend on the layout and the block's periodicity, not on which field is being iterated, so the scratch vectors this module allocates share them.

Definition at line 63 of file les.c.

{
    PetscFunctionBeginUser;
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS, PICURV_STATISTICS_MASK_FLUID, plan));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildWallModelExclusionMask()

static PetscErrorCode BuildWallModelExclusionMask ( UserCtx *  user, Vec  inclusion  )

static

Marks the near-wall cells a wall model overwrites, so the procedure skips them.

A wall model replaces the first interior cell's velocity with an analytic profile. The Germano identity evaluated there measures that imposed profile rather than resolved turbulence, so the cell is excluded from both contraction sums. This is not a modelling preference: the data at those cells is not a measurement, and averaging it in biases whatever set it belongs to. Under global averaging that set is the whole domain; under homogeneous averaging retaining the wall-normal direction it is the wall plane, whose value would otherwise be entirely imposed.

A cell left with no data receives a zero coefficient, which is the right answer where a wall model is already supplying the stress.

Only meaningful while a wall model is active; with none, nothing is imposed and every cell carries a real measurement.

Definition at line 87 of file les.c.

{
    const DMDALocalInfo info = user->info;
    const BCFace negative[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace positive[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
    const PetscInt extent[3] = {info.mx, info.my, info.mz};
    PetscReal ***mask = NULL;
    SpatialTargetPlan plan;
 
    PetscFunctionBeginUser;
    PetscCall(VecSet(inclusion, 1.0));
    PetscCall(LESResolveDomain(user, &plan));
    PetscCall(DMDAVecGetArray(user->da, inclusion, &mask));
    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++) {
        const PetscInt coord[3] = {i, j, k};
 
        for (PetscInt axis = 0; axis < 3; axis++) {
            /* Index 1 and extent-2 are the first interior cells the wall model writes. */
            if ((user->boundary_faces[negative[axis]].mathematical_type == WALL &&
                 coord[axis] == 1) ||
                (user->boundary_faces[positive[axis]].mathematical_type == WALL &&
                 coord[axis] == extent[axis] - 2)) {
                mask[k][j][i] = 0.0;
                break;
            }
        }
    }
    PetscCall(DMDAVecRestoreArray(user->da, inclusion, &mask));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LESPeriodicAxes()

static void LESPeriodicAxes ( const SimCtx *  simCtx, PetscBool  periodic[3]  )

static

Reports the block's per-axis periodicity in (xi, eta, zeta) order.

Reads the resolved global flags rather than the boundary face configuration. The two cannot disagree: the flags are derived from the boundary files, and the loader rejects a case whose opposite faces or whose blocks do not agree on periodicity. The flags are also what the DMDA itself was built from, and what the shared spatial target consults, so taking them here keeps one source of truth.

Definition at line 129 of file les.c.

{
    periodic[0] = (PetscBool)(simCtx->i_periodic != 0);
    periodic[1] = (PetscBool)(simCtx->j_periodic != 0);
    periodic[2] = (PetscBool)(simCtx->k_periodic != 0);
}

Here is the caller graph for this function:

SymTensorSelfOuter()

SymTensor SymTensorSelfOuter

(

Cmpnts

v

)

Implementation of SymTensorSelfOuter().

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
	)

Implementation of SymTensorCombine().

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

)

Implementation of SymTensorTrace().

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

)

Implementation of SymTensorDeviator().

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
	)

Implementation of SymTensorContract().

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
	)

Implementation of SymTensorTimesVector().

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

)

Implementation of SymTensorNormSq().

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
	)

Implementation of StrainRateFromGradients().

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
	)

Implementation of ComputeCellFilterWidth().

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
	)

Implementation of LeonardStress().

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
	)

Implementation of GermanoModelTensor().

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]
	)

Implementation of ResolveLESAveragingDirections().

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
	)

Implementation of ClipModelCoefficient().

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
	)

Implementation of EddyViscosityFromCoefficient().

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
	)

Implementation of WALEEddyViscosity().

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
	)

Implementation of VremanEddyViscosity().

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
	)

Implementation of SubgridKineticEnergy().

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:

LESStrainSampleIndex()

static PetscInt LESStrainSampleIndex ( PetscInt  index, PetscInt  extent, PetscBool  periodic  )

static

Maps a cell index to the index whose strain rate represents it.

The strain rate must be known one cell beyond the computed range, because the test-filter stencil reaches there. At an interior index the answer is the index itself. At the two layout-boundary planes it is not, and for opposite reasons.

On a periodic axis, index 0 and index m-1 are duplicates of the last and first interior cells. A central difference evaluated at index 0 straddles the wrap and returns nonsense, so the strain is evaluated instead at ghost index -2, which addresses the same physical cell from the interior side. Index m+1 plays the same role at the other end. These are the ghost indices the periodic cell synchronization itself copies from, and the DMDA carries the three ghost layers they need whenever any axis is periodic.

On a non-periodic axis the two planes are the solver's boundary layer, where a central difference would reach outside the domain. The nearest interior cell's strain is used instead, which is a zeroth-order extrapolation into a layer that only ever contributes to a filter average.

Definition at line 624 of file les.c.

{
    if (index == 0)          return periodic ? -2 : 1;
    if (index == extent - 1) return periodic ? extent + 1 : extent - 2;
    return index;
}

Here is the caller graph for this function:

ComputeStrainRateField()

static PetscErrorCode ComputeStrainRateField ( UserCtx *  user, Cmpnts ***  ucat, Vec  strain_diagonal, Vec  strain_offdiagonal, Vec  strain_magnitude  )

static

Fills the strain-rate scratch fields across the computed range and its halo.

Definition at line 634 of file les.c.

{
    DM                da = user->da, fda = user->fda;
    DMDALocalInfo     info;
    SpatialTargetPlan plan;
    PetscBool         periodic[3];
    Cmpnts         ***diagonal = NULL, ***offdiagonal = NULL;
    PetscReal      ***magnitude = NULL;
 
    PetscFunctionBeginUser;
 
    PetscCall(DMDAGetLocalInfo(da, &info));
    PetscCall(LESResolveDomain(user, &plan));
    LESPeriodicAxes(user->simCtx, periodic);
 
    PetscCall(DMDAVecGetArray(fda, strain_diagonal, &diagonal));
    PetscCall(DMDAVecGetArray(fda, strain_offdiagonal, &offdiagonal));
    PetscCall(DMDAVecGetArray(da, strain_magnitude, &magnitude));
 
    // One cell beyond the computed range in every direction, because that is how far
    // the 3x3x3 test-filter stencil reaches from the outermost computed cell.
    for (PetscInt k = plan.start[2] - 1; k <= plan.end[2]; k++)
    for (PetscInt j = plan.start[1] - 1; j <= plan.end[1]; j++)
    for (PetscInt i = plan.start[0] - 1; i <= plan.end[0]; i++) {
        const PetscInt si = LESStrainSampleIndex(i, info.mx, periodic[0]);
        const PetscInt sj = LESStrainSampleIndex(j, info.my, periodic[1]);
        const PetscInt sk = LESStrainSampleIndex(k, info.mz, periodic[2]);
 
        Cmpnts    dudx, dvdx, dwdx;
        SymTensor strain;
        PetscReal strain_abs;
 
        PetscCall(ComputeVectorFieldDerivatives(user, si, sj, sk, ucat, &dudx, &dvdx, &dwdx));
        PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, &strain, &strain_abs));
 
        diagonal[k][j][i].x    = strain.xx;
        diagonal[k][j][i].y    = strain.yy;
        diagonal[k][j][i].z    = strain.zz;
        offdiagonal[k][j][i].x = strain.xy;
        offdiagonal[k][j][i].y = strain.xz;
        offdiagonal[k][j][i].z = strain.yz;
        magnitude[k][j][i]     = strain_abs;
    }
 
    PetscCall(DMDAVecRestoreArray(fda, strain_diagonal, &diagonal));
    PetscCall(DMDAVecRestoreArray(fda, strain_offdiagonal, &offdiagonal));
    PetscCall(DMDAVecRestoreArray(da, strain_magnitude, &magnitude));
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LESStrainAt()

static SymTensor LESStrainAt ( Cmpnts  diagonal, Cmpnts  offdiagonal  )

inlinestatic

Rebuilds a symmetric strain tensor from its two scratch storage vectors.

Definition at line 690 of file les.c.

{
    SymTensor s;
 
    s.xx = diagonal.x;
    s.yy = diagonal.y;
    s.zz = diagonal.z;
    s.xy = offdiagonal.x;
    s.xz = offdiagonal.y;
    s.yz = offdiagonal.z;
    return s;
}

Here is the caller graph for this function:

ComputeSmagorinskyConstant()

PetscErrorCode ComputeSmagorinskyConstant

(

UserCtx *

user

)

Implementation of ComputeSmagorinskyConstant().

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

)

Implementation of ComputeEddyViscosityLES().

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

)

Implementation of LogLESDiagnostics().

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:

Variable Documentation

LES_EPSILON

const double LES_EPSILON = 1.0e-12

static

Definition at line 31 of file les.c.