test_les.c

Go to the documentation of this file.

/**
 * @file test_les.c
 * @brief C unit tests for the LES subgrid-scale closure.
 *
 * The suite is layered the way the module is. Pure kernels are checked against hand
 * computed values with no PETSc objects at all; the Germano pieces are checked on
 * fields whose answers are known in closed form; and the driver routines are checked
 * on a small DMDA fixture.
 *
 * Two of the cases exist because of specific defects and are worth naming. The
 * "filtered-product" case pins the distinction between the filter of a product and
 * the product of filtered factors, which the model tensor once collapsed. The
 * "duplicate-plane" case pins the requirement that a spatial average count each
 * physical cell once on a periodic block, where the two layout boundary planes are
 * copies of interior cells.
 */
 
#include "test_support.h"
 
#include "les.h"
#include "setup.h"
#include "statistics_target.h"
 
/*================================================================================*
 *                              KERNEL-LEVEL TESTS                                *
 *================================================================================*/
 
/** @brief Tests the symmetric-tensor primitives against hand-computed values. */
static PetscErrorCode TestSymTensorAlgebra(void)
{
    const Cmpnts    v = {1.0, 2.0, 3.0};
    const SymTensor a = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
    const SymTensor b = {2.0, 0.0, 0.0, 2.0, 0.0, 2.0};
    SymTensor       outer, deviator, combined;
 
    PetscFunctionBeginUser;
 
    outer = SymTensorSelfOuter(v);
    PetscCall(PicurvAssertRealNear(1.0, outer.xx, 1.0e-14, "self outer product should square the x component"));
    PetscCall(PicurvAssertRealNear(2.0, outer.xy, 1.0e-14, "self outer product should cross x and y"));
    PetscCall(PicurvAssertRealNear(6.0, outer.yz, 1.0e-14, "self outer product should cross y and z"));
    PetscCall(PicurvAssertRealNear(9.0, outer.zz, 1.0e-14, "self outer product should square the z component"));
 
    /* trace(a) = 1 + 4 + 6 = 11 */
    PetscCall(PicurvAssertRealNear(11.0, SymTensorTrace(a), 1.0e-14, "trace should sum the diagonal"));
 
    deviator = SymTensorDeviator(a);
    PetscCall(PicurvAssertRealNear(0.0, SymTensorTrace(deviator), 1.0e-13,
                                   "the deviator should be trace free"));
    PetscCall(PicurvAssertRealNear(a.xy, deviator.xy, 1.0e-14,
                                   "the deviator should leave off-diagonal components alone"));
 
    /* a:b doubles the off-diagonal terms; b is diagonal, so a:b = 2*(1 + 4 + 6) = 22. */
    PetscCall(PicurvAssertRealNear(22.0, SymTensorContract(a, b), 1.0e-13,
                                   "contraction against a diagonal tensor should scale its trace"));
    /* a:a = 1 + 16 + 36 + 2*(4 + 9 + 25) = 53 + 76 = 129 */
    PetscCall(PicurvAssertRealNear(129.0, SymTensorNormSq(a), 1.0e-13,
                                   "the squared norm must count off-diagonal components twice"));
 
    combined = SymTensorCombine(2.0, a, -1.0, b);
    PetscCall(PicurvAssertRealNear(0.0, combined.xx, 1.0e-14, "combine should scale and subtract"));
    PetscCall(PicurvAssertRealNear(4.0, combined.xy, 1.0e-14, "combine should scale off-diagonal terms"));
 
    PetscFunctionReturn(0);
}
 
/** @brief Tests strain-rate assembly and its magnitude for a known velocity gradient. */
static PetscErrorCode TestStrainRateFromGradients(void)
{
    /* du/dy = 2 and every other derivative zero: pure shear. */
    const Cmpnts dudx = {0.0, 2.0, 0.0};
    const Cmpnts dvdx = {0.0, 0.0, 0.0};
    const Cmpnts dwdx = {0.0, 0.0, 0.0};
    SymTensor    strain;
    PetscReal    magnitude;
 
    PetscFunctionBeginUser;
    PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, &strain, &magnitude));
 
    PetscCall(PicurvAssertRealNear(1.0, strain.xy, 1.0e-14,
                                   "pure shear should give S_xy = half the velocity gradient"));
    PetscCall(PicurvAssertRealNear(0.0, strain.xx, 1.0e-14, "pure shear should have no normal strain"));
    PetscCall(PicurvAssertRealNear(0.0, SymTensorTrace(strain), 1.0e-14,
                                   "an incompressible gradient should give a trace-free strain"));
    /* |S| = sqrt(2 S_ij S_ij) = sqrt(2 * 2 * 1) = 2 */
    PetscCall(PicurvAssertRealNear(2.0, magnitude, 1.0e-13,
                                   "the strain magnitude must count both off-diagonal entries"));
 
    PetscFunctionReturn(0);
}
 
/** @brief Tests that the three filter-width models separate on an anisotropic cell. */
static PetscErrorCode TestFilterWidthModelsSeparateOnStretchedCell(void)
{
    /* Cartesian-aligned metrics for a cell of size 4 x 1 x 1: the face area vectors
       are the products of the two spanning edges, and the Jacobian is 1/volume. */
    const Cmpnts    csi = {1.0, 0.0, 0.0};
    const Cmpnts    eta = {0.0, 4.0, 0.0};
    const Cmpnts    zet = {0.0, 0.0, 4.0};
    const PetscReal aj  = 0.25;
    PetscReal       cube_root = 0.0, geometric = 0.0, max_edge = 0.0;
 
    PetscFunctionBeginUser;
 
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_CUBE_ROOT_VOLUME, aj, csi, eta, zet, &cube_root));
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_GEOMETRIC_MEAN, aj, csi, eta, zet, &geometric));
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_MAX_EDGE, aj, csi, eta, zet, &max_edge));
 
    /* Volume is 4, so the cube-root model gives 4^(1/3). */
    PetscCall(PicurvAssertRealNear(PetscPowReal(4.0, 1.0 / 3.0), cube_root, 1.0e-12,
                                   "the cube-root model should depend only on cell volume"));
    PetscCall(PicurvAssertRealNear(PetscPowReal(4.0, 1.0 / 3.0), geometric, 1.0e-12,
                                   "the geometric mean of the extents should recover the same volume"));
    PetscCall(PicurvAssertRealNear(4.0, max_edge, 1.0e-12,
                                   "the max-edge model should report the long direction"));
    PetscCall(PicurvAssertBool((PetscBool)(max_edge > cube_root),
                               "a stretched cell should make the max-edge width the larger of the two"));
 
    PetscFunctionReturn(0);
}
 
/** @brief Face-area vectors and inverse volume of the parallelepiped spanned by a, b, c. */
static void ParallelepipedMetrics(Cmpnts a, Cmpnts b, Cmpnts c,
                                  Cmpnts *csi, Cmpnts *eta, Cmpnts *zet, PetscReal *aj)
{
    /* The same orientation convention as ComputeFaceMetrics(): each face-area vector
       is the cross product of the two edges spanning that face. */
    csi->x = b.y*c.z - b.z*c.y; csi->y = b.z*c.x - b.x*c.z; csi->z = b.x*c.y - b.y*c.x;
    eta->x = c.y*a.z - c.z*a.y; eta->y = c.z*a.x - c.x*a.z; eta->z = c.x*a.y - c.y*a.x;
    zet->x = a.y*b.z - a.z*b.y; zet->y = a.z*b.x - a.x*b.z; zet->z = a.x*b.y - a.y*b.x;
    *aj = 1.0 / (a.x*csi->x + a.y*csi->y + a.z*csi->z);
}
 
/**
 * @brief Tests that a filter width is a property of the cell, not of its orientation.
 *
 * @details The widths used to come from the Cartesian components of the cell diagonal.
 *          On the Humphrey bend an identical wall cell got a geometric-mean width three
 *          times larger, and so nine times the eddy viscosity, merely because the bend
 *          had turned it through 45 degrees relative to the x and y axes.
 */
static PetscErrorCode TestFilterWidthIsIndependentOfCellOrientation(void)
{
    /* A thin, long cell - 0.25 x 1 x 8 - first axis-aligned, then turned about z by
       45 degrees and tilted about x by 30 degrees. */
    const PetscReal c45 = PetscCosReal(PETSC_PI / 4.0), s45 = PetscSinReal(PETSC_PI / 4.0);
    const PetscReal c30 = PetscCosReal(PETSC_PI / 6.0), s30 = PetscSinReal(PETSC_PI / 6.0);
    const Cmpnts    aligned[3] = {{0.25, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 8.0}};
    const LESFilterWidthModel models[3] = {
        LES_FILTER_WIDTH_CUBE_ROOT_VOLUME, LES_FILTER_WIDTH_GEOMETRIC_MEAN, LES_FILTER_WIDTH_MAX_EDGE};
    Cmpnts          rotated[3];
    Cmpnts          csi, eta, zet;
    PetscReal       aj, reference[3], turned[3];
 
    PetscFunctionBeginUser;
    for (PetscInt e = 0; e < 3; ++e) {
        /* Rz(45) then Rx(30). */
        const Cmpnts v = aligned[e];
        const Cmpnts z = {c45*v.x - s45*v.y, s45*v.x + c45*v.y, v.z};
        rotated[e] = (Cmpnts){z.x, c30*z.y - s30*z.z, s30*z.y + c30*z.z};
    }
 
    ParallelepipedMetrics(aligned[0], aligned[1], aligned[2], &csi, &eta, &zet, &aj);
    for (PetscInt m = 0; m < 3; ++m) PetscCall(ComputeCellFilterWidth(models[m], aj, csi, eta, zet, &reference[m]));
    ParallelepipedMetrics(rotated[0], rotated[1], rotated[2], &csi, &eta, &zet, &aj);
    for (PetscInt m = 0; m < 3; ++m) PetscCall(ComputeCellFilterWidth(models[m], aj, csi, eta, zet, &turned[m]));
 
    PetscCall(PicurvAssertRealNear(reference[0], turned[0], 1.0e-12, "cube-root width is orientation independent"));
    PetscCall(PicurvAssertRealNear(reference[1], turned[1], 1.0e-12, "geometric-mean width is orientation independent"));
    PetscCall(PicurvAssertRealNear(reference[2], turned[2], 1.0e-12, "max-edge width is orientation independent"));
    PetscCall(PicurvAssertRealNear(8.0, turned[2], 1.0e-12, "max-edge width is the cell's longest extent"));
    PetscCall(PicurvAssertRealNear(PetscCbrtReal(2.0), turned[1], 1.0e-12,
                                   "geometric-mean width of a 0.25 x 1 x 8 cell is the cube root of its volume"));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests the Scotti aspect-ratio correction against its closed form.
 *
 * @details Scotti, Meneveau & Lilly (1993) multiply the cube root of the extents by
 *          f(a1, a2) = cosh(sqrt((4/27)(ln^2 a1 - ln a1 ln a2 + ln^2 a2))), with a1, a2
 *          the two shorter extents over the longest. f is exactly 1 for a cube, so the
 *          width must equal the cube root there, and it must be unchanged by rotating
 *          the cell.
 */
static PetscErrorCode TestScottiFilterWidthMatchesClosedForm(void)
{
    const Cmpnts cube[3]  = {{1.5, 0.0, 0.0}, {0.0, 1.5, 0.0}, {0.0, 0.0, 1.5}};
    const Cmpnts brick[3] = {{1.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 8.0}};
    const PetscReal c45 = PetscCosReal(PETSC_PI / 4.0), s45 = PetscSinReal(PETSC_PI / 4.0);
    Cmpnts    turned[3], csi, eta, zet;
    PetscReal aj, delta, delta_turned, cube_root;
 
    PetscFunctionBeginUser;
    ParallelepipedMetrics(cube[0], cube[1], cube[2], &csi, &eta, &zet, &aj);
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_SCOTTI, aj, csi, eta, zet, &delta));
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_CUBE_ROOT_VOLUME, aj, csi, eta, zet, &cube_root));
    PetscCall(PicurvAssertRealNear(cube_root, delta, 1.0e-12, "the Scotti correction is 1 for a cube"));
 
    /* 1 x 2 x 8: a1 = 1/8, a2 = 1/4, f = 1.2596408907321084. */
    ParallelepipedMetrics(brick[0], brick[1], brick[2], &csi, &eta, &zet, &aj);
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_SCOTTI, aj, csi, eta, zet, &delta));
    PetscCall(PicurvAssertRealNear(3.1740961470834224, delta, 1.0e-12,
                                   "Scotti width of a 1 x 2 x 8 cell matches the closed form"));
 
    for (PetscInt e = 0; e < 3; ++e) {
        const Cmpnts v = brick[e];
        turned[e] = (Cmpnts){c45 * v.x - s45 * v.z, v.y, s45 * v.x + c45 * v.z};
    }
    ParallelepipedMetrics(turned[0], turned[1], turned[2], &csi, &eta, &zet, &aj);
    PetscCall(ComputeCellFilterWidth(LES_FILTER_WIDTH_SCOTTI, aj, csi, eta, zet, &delta_turned));
    PetscCall(PicurvAssertRealNear(delta, delta_turned, 1.0e-12, "the Scotti width is orientation independent"));
    PetscFunctionReturn(0);
}
 
/** @brief Rotates a vector by Rz(40 deg) then Rx(25 deg). */
static Cmpnts RotateForInvarianceTest(Cmpnts v)
{
    const PetscReal cz = PetscCosReal(40.0 * PETSC_PI / 180.0), sz = PetscSinReal(40.0 * PETSC_PI / 180.0);
    const PetscReal cx = PetscCosReal(25.0 * PETSC_PI / 180.0), sx = PetscSinReal(25.0 * PETSC_PI / 180.0);
    const Cmpnts z = {cz * v.x - sz * v.y, sz * v.x + cz * v.y, v.z};
    return (Cmpnts){z.x, cx * z.y - sx * z.z, sx * z.y + cx * z.z};
}
 
/** @brief Rotates a velocity gradient g_ij = du_i/dx_j as R g R^T, rows passed as dudx, dvdx, dwdx. */
static void RotateGradientForInvarianceTest(const Cmpnts rows[3], Cmpnts out[3])
{
    Cmpnts column_rotated[3], cols[3];
    /* g R^T: rotate each row as a covector, i.e. rotate every row vector. */
    for (PetscInt a = 0; a < 3; ++a) column_rotated[a] = RotateForInvarianceTest(rows[a]);
    /* R (g R^T): rotate the columns, which are the components across rows. */
    for (PetscInt b = 0; b < 3; ++b) {
        const Cmpnts column = {b == 0 ? column_rotated[0].x : (b == 1 ? column_rotated[0].y : column_rotated[0].z),
                               b == 0 ? column_rotated[1].x : (b == 1 ? column_rotated[1].y : column_rotated[1].z),
                               b == 0 ? column_rotated[2].x : (b == 1 ? column_rotated[2].y : column_rotated[2].z)};
        cols[b] = RotateForInvarianceTest(column);
    }
    for (PetscInt a = 0; a < 3; ++a)
        out[a] = (Cmpnts){a == 0 ? cols[0].x : (a == 1 ? cols[0].y : cols[0].z),
                          a == 0 ? cols[1].x : (a == 1 ? cols[1].y : cols[1].z),
                          a == 0 ? cols[2].x : (a == 1 ? cols[2].y : cols[2].z)};
}
 
/**
 * @brief Tests the WALE and Vreman kernels against independent evaluations.
 *
 * @details Reference values come from a separate numpy evaluation of the published
 *          formulas for the gradient below. Both models must vanish in pure shear, which
 *          is the property that lets them stay silent in a laminar shear layer, and both
 *          must be unchanged when the velocity gradient and the cell rotate together.
 */
static PetscErrorCode TestWALEAndVremanKernels(void)
{
    const Cmpnts g[3] = {{0.3, 1.2, -0.4}, {0.7, -0.5, 0.9}, {-0.2, 0.6, 0.2}};
    const Cmpnts shear[3] = {{0.0, 2.5, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
    const Cmpnts aligned_edges[3] = {{0.1, 0.0, 0.0}, {0.0, 0.2, 0.0}, {0.0, 0.0, 0.8}};
    Cmpnts rotated_edges[3], rotated_g[3], csi, eta, zet, recovered[3];
    PetscReal aj;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvAssertRealNear(0.0002799988905128629, WALEEddyViscosity(g[0], g[1], g[2], 0.1, 0.5),
                                   1.0e-15, "WALE matches the independent evaluation"));
    PetscCall(PicurvAssertRealNear(0.007234428329928329,
                                   VremanEddyViscosity(g[0], g[1], g[2], aligned_edges, 0.07),
                                   1.0e-15, "Vreman on aligned edges matches the Cartesian formula"));
    PetscCall(PicurvAssertRealNear(0.0, WALEEddyViscosity(shear[0], shear[1], shear[2], 0.1, 0.5),
                                   1.0e-15, "WALE vanishes in pure shear"));
    PetscCall(PicurvAssertRealNear(0.0, VremanEddyViscosity(shear[0], shear[1], shear[2], aligned_edges, 0.07),
                                   1.0e-15, "Vreman vanishes in pure shear"));
 
    for (PetscInt m = 0; m < 3; ++m) rotated_edges[m] = RotateForInvarianceTest(aligned_edges[m]);
    RotateGradientForInvarianceTest(g, rotated_g);
    PetscCall(PicurvAssertRealNear(VremanEddyViscosity(g[0], g[1], g[2], aligned_edges, 0.07),
                                   VremanEddyViscosity(rotated_g[0], rotated_g[1], rotated_g[2], rotated_edges, 0.07),
                                   1.0e-14, "Vreman is unchanged when the flow and the cell rotate together"));
    PetscCall(PicurvAssertRealNear(WALEEddyViscosity(g[0], g[1], g[2], 0.1, 0.5),
                                   WALEEddyViscosity(rotated_g[0], rotated_g[1], rotated_g[2], 0.1, 0.5),
                                   1.0e-15, "WALE is unchanged when the flow rotates"));
 
    /* The edge vectors recovered from a turned cell's metrics are its actual edges. */
    ParallelepipedMetrics(rotated_edges[0], rotated_edges[1], rotated_edges[2], &csi, &eta, &zet, &aj);
    PetscCall(ComputeCellEdgeVectors(aj, csi, eta, zet, recovered));
    for (PetscInt m = 0; m < 3; ++m) {
        PetscCall(PicurvAssertRealNear(rotated_edges[m].x, recovered[m].x, 1.0e-14, "edge x recovered from metrics"));
        PetscCall(PicurvAssertRealNear(rotated_edges[m].y, recovered[m].y, 1.0e-14, "edge y recovered from metrics"));
        PetscCall(PicurvAssertRealNear(rotated_edges[m].z, recovered[m].z, 1.0e-14, "edge z recovered from metrics"));
    }
    PetscFunctionReturn(0);
}
 
/** @brief Tests that the Leonard stress vanishes on a uniform velocity field. */
static PetscErrorCode TestLeonardStressVanishesOnUniformFlow(void)
{
    const Cmpnts    velocity = {1.5, -0.5, 2.0};
    const SymTensor product  = SymTensorSelfOuter(velocity);
    SymTensor       leonard;
 
    PetscFunctionBeginUser;
    /* Filtering a constant returns the constant, so both filtered inputs are exact. */
    leonard = LeonardStress(velocity, product);
 
    PetscCall(PicurvAssertRealNear(0.0, SymTensorNormSq(leonard), 1.0e-24,
                                   "a uniform field carries no stress between the two filter widths"));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests the model tensor where the two filter terms provably coincide.
 *
 * On a constant strain field the test filter is the identity, so the filter of the
 * product equals the product of the filtered factors and the model tensor collapses
 * to a closed form, `-2 Delta^2 (alpha - 1) |S| S_ij`, made trace free. This is the
 * one configuration in which the corrected tensor and the collapsed expression that
 * preceded it agree, which is what makes it a clean analytic anchor.
 */
static PetscErrorCode TestGermanoModelTensorOnConstantStrain(void)
{
    const PetscReal delta = 0.5;
    const PetscReal alpha = 4.0;
    Cmpnts          dudx = {0.0, 2.0, 0.0}, dvdx = {0.0, 0.0, 0.0}, dwdx = {0.0, 0.0, 0.0};
    SymTensor       strain, product, model, expected;
    PetscReal       magnitude;
 
    PetscFunctionBeginUser;
    PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, &strain, &magnitude));
 
    product = SymTensorCombine(magnitude, strain, 0.0, strain);
    model   = GermanoModelTensor(delta, alpha, magnitude, strain, product);
 
    expected = SymTensorDeviator(SymTensorCombine(-2.0 * delta * delta * (alpha - 1.0) * magnitude,
                                                  strain, 0.0, strain));
 
    PetscCall(PicurvAssertRealNear(expected.xy, model.xy, 1.0e-13,
                                   "constant strain should give the closed-form model tensor"));
    PetscCall(PicurvAssertRealNear(0.0, SymTensorTrace(model), 1.0e-13,
                                   "the model tensor must be trace free"));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the model tensor uses the filtered product, not the filtered factors.
 *
 * Builds a stencil whose strain magnitude varies across it, so the test filter of
 * `|S| S_ij` genuinely differs from `|S|^ S^_ij`. The model tensor must be built from
 * the former. Collapsing it onto the latter yields `-2 Delta^2 (alpha - 1) |S|^ S^_ij`,
 * which is well scaled, responds to the flow, and is not the Germano-Lilly tensor;
 * this case fails if the two are ever conflated again.
 */
static PetscErrorCode TestGermanoModelTensorUsesFilteredProduct(void)
{
    const LESTestFilterKernel kernel = LES_TEST_FILTER_VOLUME_WEIGHTED_BOX;
    const PetscReal           delta = 1.0;
    const PetscReal           alpha = 4.0;
    double                    magnitude_stencil[3][3][3], weights[3][3][3];
    SymTensor                 strain_stencil[3][3][3], product_stencil[3][3][3];
    SymTensor                 strain_filtered, product_filtered, model, collapsed;
    PetscReal                 magnitude_filtered, separation;
 
    PetscFunctionBeginUser;
 
    for (PetscInt k = 0; k < 3; ++k)
    for (PetscInt j = 0; j < 3; ++j)
    for (PetscInt i = 0; i < 3; ++i) {
        /* A strain that grows across the stencil: filtering no longer commutes with
           forming the product, which is exactly the regime the identity relies on. */
        const PetscReal shear = 1.0 + 0.5 * (PetscReal)(i + j + k);
        Cmpnts          dudx = {0.0, 2.0 * shear, 0.0};
        Cmpnts          dvdx = {0.0, 0.0, 0.0};
        Cmpnts          dwdx = {0.0, 0.0, 0.0};
        SymTensor       strain;
        PetscReal       magnitude;
 
        PetscCall(StrainRateFromGradients(dudx, dvdx, dwdx, &strain, &magnitude));
        strain_stencil[k][j][i]    = strain;
        magnitude_stencil[k][j][i] = magnitude;
        product_stencil[k][j][i]   = SymTensorCombine(magnitude, strain, 0.0, strain);
        weights[k][j][i]           = 1.0;
    }
 
    magnitude_filtered = ApplyLESTestFilter(kernel, magnitude_stencil, weights);
    PetscCall(ApplyLESTestFilterSymTensor(kernel, strain_stencil, weights, &strain_filtered));
    PetscCall(ApplyLESTestFilterSymTensor(kernel, product_stencil, weights, &product_filtered));
 
    /* The filter of the product must not equal the product of the filtered factors. */
    separation = PetscAbsReal(product_filtered.xy - magnitude_filtered * strain_filtered.xy);
    PetscCall(PicurvAssertBool((PetscBool)(separation > 1.0e-6),
                               "a varying strain field must separate the filtered product from "
                               "the product of the filtered factors"));
 
    model     = GermanoModelTensor(delta, alpha, magnitude_filtered, strain_filtered, product_filtered);
    collapsed = SymTensorDeviator(SymTensorCombine(-2.0 * delta * delta * (alpha - 1.0) * magnitude_filtered,
                                                   strain_filtered, 0.0, strain_filtered));
 
    PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(model.xy - collapsed.xy) > 1.0e-6),
                               "the model tensor must be built from the filtered product, not from "
                               "the product of separately filtered quantities"));
    PetscFunctionReturn(0);
}
 
/** @brief Tests each limiting mode, including the sign that carries backscatter. */
static PetscErrorCode TestClipModelCoefficientModes(void)
{
    LESConfig config;
    PetscBool limited = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(LESConfigSetDefaults(&config));
    config.max_cs = 0.2; /* ceiling on Cs, so the coefficient ceiling is 0.04 */
 
    config.clip_mode = LES_CLIP_CLAMP;
    PetscCall(PicurvAssertRealNear(0.0, ClipModelCoefficient(-0.01, &config, &limited), 1.0e-14,
                                   "clamping should discard a negative coefficient"));
    PetscCall(PicurvAssertBool(limited, "clamping a negative coefficient counts as limiting"));
    PetscCall(PicurvAssertRealNear(0.04, ClipModelCoefficient(1.0, &config, &limited), 1.0e-14,
                                   "clamping should cap the coefficient at max_cs squared"));
    PetscCall(PicurvAssertRealNear(0.01, ClipModelCoefficient(0.01, &config, &limited), 1.0e-14,
                                   "clamping should pass an admissible coefficient through"));
    PetscCall(PicurvAssertBool((PetscBool)(!limited), "an admissible coefficient is not limited"));
 
    config.clip_mode = LES_CLIP_CLIP_NEGATIVE;
    PetscCall(PicurvAssertRealNear(0.0, ClipModelCoefficient(-0.01, &config, &limited), 1.0e-14,
                                   "clip_negative should discard a negative coefficient"));
    PetscCall(PicurvAssertRealNear(1.0, ClipModelCoefficient(1.0, &config, &limited), 1.0e-14,
                                   "clip_negative should impose no ceiling"));
 
    config.clip_mode = LES_CLIP_NONE;
    PetscCall(PicurvAssertRealNear(-0.01, ClipModelCoefficient(-0.01, &config, &limited), 1.0e-14,
                                   "no limiting should preserve the sign that carries backscatter"));
    PetscCall(PicurvAssertBool((PetscBool)(!limited), "an unlimited coefficient is not modified"));
    PetscFunctionReturn(0);
}
 
/** @brief Tests eddy-viscosity assembly and the total-viscosity floor. */
static PetscErrorCode TestEddyViscosityFloorBoundsTotalViscosity(void)
{
    const PetscReal molecular = 0.01;
 
    PetscFunctionBeginUser;
 
    /* C = 0.04, Delta = 2, |S| = 3  ->  nu_t = 0.04 * 4 * 3 = 0.48 */
    PetscCall(PicurvAssertRealNear(0.48, EddyViscosityFromCoefficient(0.04, 2.0, 3.0, molecular, 0.0),
                                   1.0e-13, "eddy viscosity should be the coefficient times Delta^2 |S|"));
 
    /* A backscattering coefficient may drive nu_t negative, but only down to the point
       where the total viscosity would vanish. */
    PetscCall(PicurvAssertRealNear(-molecular,
                                   EddyViscosityFromCoefficient(-100.0, 2.0, 3.0, molecular, 0.0),
                                   1.0e-14, "the floor should stop nu_t at minus the molecular value"));
    PetscCall(PicurvAssertRealNear(-0.5 * molecular,
                                   EddyViscosityFromCoefficient(-100.0, 2.0, 3.0, molecular, 0.5),
                                   1.0e-14, "a ratio of one half should keep half the molecular viscosity"));
    PetscCall(PicurvAssertRealNear(-0.006,
                                   EddyViscosityFromCoefficient(-0.0005, 2.0, 3.0, molecular, 0.0),
                                   1.0e-14, "a mild backscatter should survive the floor untouched"));
    PetscFunctionReturn(0);
}
 
/** @brief Tests the Yoshizawa subgrid kinetic energy against its closed form. */
static PetscErrorCode TestSubgridKineticEnergy(void)
{
    PetscFunctionBeginUser;
    /* 2 * 0.09 * 4 * 9 = 6.48 */
    PetscCall(PicurvAssertRealNear(6.48, SubgridKineticEnergy(0.09, 2.0, 3.0), 1.0e-13,
                                   "the Yoshizawa relation should scale with Delta^2 |S|^2"));
    PetscFunctionReturn(0);
}
 
/*================================================================================*
 *                            FIXTURE-LEVEL TESTS                                 *
 *================================================================================*/
 
/**
 * @brief Fills a ghosted local scalar field over the whole local array.
 *
 * Writes the halo as well as the owned cells, so a test can poison the layout
 * boundary planes and observe whether an averaging routine counts them.
 */
static PetscErrorCode FillLocalScalar(UserCtx *user, Vec local, PetscReal value)
{
    PetscFunctionBeginUser;
    PetscCall(VecSet(local, value));
    (void)user;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Writes a scalar into every cell the closure treats as owned and interior.
 */
static PetscErrorCode SetInteriorScalar(UserCtx *user, Vec local, PetscReal value)
{
    DMDALocalInfo info;
    PetscReal  ***array = NULL;
    PetscInt      xs, xe, ys, ye, zs, ze;
 
    PetscFunctionBeginUser;
    PetscCall(DMDAGetLocalInfo(user->da, &info));
    xs = (info.xs == 0) ? 1 : info.xs;
    ys = (info.ys == 0) ? 1 : info.ys;
    zs = (info.zs == 0) ? 1 : info.zs;
    xe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
    ye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
    ze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
 
    PetscCall(DMDAVecGetArray(user->da, local, &array));
    for (PetscInt k = zs; k < ze; ++k)
    for (PetscInt j = ys; j < ye; ++j)
    for (PetscInt i = xs; i < xe; ++i) array[k][j][i] = value;
    PetscCall(DMDAVecRestoreArray(user->da, local, &array));
    PetscFunctionReturn(0);
}
 
/** @brief Reads one cell of a ghosted local scalar field. */
static PetscErrorCode ReadLocalScalar(UserCtx *user, Vec local, PetscInt i, PetscInt j, PetscInt k,
                                      PetscReal *value)
{
    PetscReal ***array = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->da, local, &array));
    *value = array[k][j][i];
    PetscCall(DMDAVecRestoreArrayRead(user->da, local, &array));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Declares periodic boundary pairs on the axes a test wants homogeneous.
 *
 * The shared fixture builds a periodic DMDA from the `SimCtx` flags but leaves the
 * boundary face configuration untouched. The periodic field synchronization reads the
 * faces rather than the flags, so a test that synchronizes has to declare them. The
 * closure's own periodicity questions are answered from the flags, so tests that only
 * ask those do not call this.
 */
static void DeclarePeriodicFaces(UserCtx *user, PetscBool xi, PetscBool eta, PetscBool zeta)
{
    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 PetscBool selected[3] = {xi, eta, zeta};
 
    for (PetscInt axis = 0; axis < 3; ++axis) {
        if (!selected[axis]) continue;
        user->boundary_faces[negative[axis]].mathematical_type = PERIODIC;
        user->boundary_faces[positive[axis]].mathematical_type = PERIODIC;
    }
}
 
/** @brief Tests that an empty direction set reproduces the pointwise local model. */
static PetscErrorCode TestAverageRatioLocalIsPointwise(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Vec          numerator, denominator, ratio;
    SpatialTargetPlan plan;
    PetscReal ***num = NULL;
    PetscReal    value = 0.0;
    const PetscBool none[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->lNvert, 0.0));
    PetscCall(VecSet(user->lAj, 1.0));
 
    PetscCall(VecDuplicate(user->lNvert, &numerator));
    PetscCall(VecDuplicate(user->lNvert, &denominator));
    PetscCall(VecDuplicate(user->lNvert, &ratio));
    PetscCall(FillLocalScalar(user, denominator, 4.0));
    PetscCall(FillLocalScalar(user, ratio, 0.0));
    PetscCall(FillLocalScalar(user, numerator, 0.0));
 
    /* A numerator that varies cell by cell, so a pointwise result cannot be confused
       with an averaged one. */
    PetscCall(DMDAVecGetArray(user->da, numerator, &num));
    for (PetscInt k = 1; k < 6; ++k)
    for (PetscInt j = 1; j < 6; ++j)
    for (PetscInt i = 1; i < 6; ++i) num[k][j][i] = (PetscReal)(i + j + k);
    PetscCall(DMDAVecRestoreArray(user->da, numerator, &num));
 
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvSpatialRatioAverage(user, &plan, numerator, denominator, NULL, none,
                                        PETSC_COMM_WORLD, ratio, NULL));
 
    PetscCall(ReadLocalScalar(user, ratio, 2, 3, 4, &value));
    PetscCall(PicurvAssertRealNear((2.0 + 3.0 + 4.0) / 4.0, value, 1.0e-13,
                                   "with no averaging direction the quotient should stay pointwise"));
 
    PetscCall(VecDestroy(&numerator));
    PetscCall(VecDestroy(&denominator));
    PetscCall(VecDestroy(&ratio));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that averaging divides summed numerators by summed denominators.
 *
 * The distinction matters: the mean of the quotients differs from the quotient of the
 * means, and Lilly's least-squares closure calls for the latter. Separating them needs
 * a denominator that varies too. Three cells carry 1/1 and two carry 6/2, so every
 * pointwise quotient is 1 or 3 and their mean is 9/5, while the ratio of the sums is
 * 15/7. Only the second answer can come from averaging the two fields first.
 */
static PetscErrorCode TestAverageRatioDividesSummedFields(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Vec          numerator, denominator, ratio;
    SpatialTargetPlan plan;
    PetscReal ***num = NULL, ***den = NULL;
    PetscReal    value = 0.0;
    const PetscBool all[3] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->lNvert, 0.0));
    PetscCall(VecSet(user->lAj, 1.0));
 
    PetscCall(VecDuplicate(user->lNvert, &numerator));
    PetscCall(VecDuplicate(user->lNvert, &denominator));
    PetscCall(VecDuplicate(user->lNvert, &ratio));
    PetscCall(FillLocalScalar(user, numerator, 0.0));
    PetscCall(FillLocalScalar(user, denominator, 0.0));
    PetscCall(FillLocalScalar(user, ratio, 0.0));
 
    PetscCall(DMDAVecGetArray(user->da, numerator, &num));
    PetscCall(DMDAVecGetArray(user->da, denominator, &den));
    for (PetscInt k = 1; k < 6; ++k)
    for (PetscInt j = 1; j < 6; ++j)
    for (PetscInt i = 1; i < 6; ++i) {
        const PetscBool low = (PetscBool)(i < 4);
 
        num[k][j][i] = low ? 1.0 : 6.0;
        den[k][j][i] = low ? 1.0 : 2.0;
    }
    PetscCall(DMDAVecRestoreArray(user->da, numerator, &num));
    PetscCall(DMDAVecRestoreArray(user->da, denominator, &den));
 
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvSpatialRatioAverage(user, &plan, numerator, denominator, NULL, all,
                                        PETSC_COMM_WORLD, ratio, NULL));
 
    /* 15/7, the ratio of the sums. Averaging the pointwise quotients would give 9/5. */
    PetscCall(ReadLocalScalar(user, ratio, 2, 3, 4, &value));
    PetscCall(PicurvAssertRealNear(15.0 / 7.0, value, 1.0e-13,
                                   "the global average should divide summed numerators by summed "
                                   "denominators, not average the pointwise quotients"));
    PetscCall(ReadLocalScalar(user, ratio, 4, 2, 1, &value));
    PetscCall(PicurvAssertRealNear(15.0 / 7.0, value, 1.0e-13,
                                   "a global average should be the same number in every cell"));
 
    PetscCall(VecDestroy(&numerator));
    PetscCall(VecDestroy(&denominator));
    PetscCall(VecDestroy(&ratio));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that averaging over two directions leaves a profile along the third.
 */
static PetscErrorCode TestAverageRatioRetainsUnaveragedDirection(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Vec          numerator, denominator, ratio;
    SpatialTargetPlan plan;
    PetscReal ***num = NULL;
    PetscReal    at_low = 0.0, at_high = 0.0;
    /* Average over xi and zeta, retain eta: the channel-flow arrangement. */
    const PetscBool ik[3] = {PETSC_TRUE, PETSC_FALSE, PETSC_TRUE};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->lNvert, 0.0));
    PetscCall(VecSet(user->lAj, 1.0));
 
    PetscCall(VecDuplicate(user->lNvert, &numerator));
    PetscCall(VecDuplicate(user->lNvert, &denominator));
    PetscCall(VecDuplicate(user->lNvert, &ratio));
    PetscCall(FillLocalScalar(user, numerator, 0.0));
    PetscCall(FillLocalScalar(user, denominator, 0.0));
    PetscCall(FillLocalScalar(user, ratio, 0.0));
    PetscCall(SetInteriorScalar(user, denominator, 1.0));
 
    /* Depends on eta alone, so averaging over xi and zeta must leave it untouched. */
    PetscCall(DMDAVecGetArray(user->da, numerator, &num));
    for (PetscInt k = 1; k < 6; ++k)
    for (PetscInt j = 1; j < 6; ++j)
    for (PetscInt i = 1; i < 6; ++i) num[k][j][i] = (PetscReal)j;
    PetscCall(DMDAVecRestoreArray(user->da, numerator, &num));
 
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvSpatialRatioAverage(user, &plan, numerator, denominator, NULL, ik,
                                        PETSC_COMM_WORLD, ratio, NULL));
 
    PetscCall(ReadLocalScalar(user, ratio, 2, 1, 3, &at_low));
    PetscCall(ReadLocalScalar(user, ratio, 4, 5, 2, &at_high));
    PetscCall(PicurvAssertRealNear(1.0, at_low, 1.0e-13,
                                   "averaging over xi and zeta should preserve the eta profile"));
    PetscCall(PicurvAssertRealNear(5.0, at_high, 1.0e-13,
                                   "each retained eta plane should keep its own average"));
 
    PetscCall(VecDestroy(&numerator));
    PetscCall(VecDestroy(&denominator));
    PetscCall(VecDestroy(&ratio));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that a spatial average ignores the periodic duplicate planes.
 *
 * On a periodic block, index 0 and index `m-1` hold copies of interior cells. Counting
 * them would weight those cells twice. The test poisons both planes and requires the
 * average not to move.
 */
static PetscErrorCode TestAverageRatioIgnoresPeriodicDuplicatePlanes(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Vec          numerator, denominator, ratio;
    SpatialTargetPlan plan;
    PetscReal ***num = NULL;
    DMDALocalInfo info;
    PetscReal    clean = 0.0, poisoned = 0.0;
    const PetscBool all[3] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 6, 6, 6,
                                                         PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    DeclarePeriodicFaces(user, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE);
    PetscCall(DMDAGetLocalInfo(user->da, &info));
    PetscCall(VecSet(user->lNvert, 0.0));
    PetscCall(VecSet(user->lAj, 1.0));
 
    PetscCall(VecDuplicate(user->lNvert, &numerator));
    PetscCall(VecDuplicate(user->lNvert, &denominator));
    PetscCall(VecDuplicate(user->lNvert, &ratio));
    PetscCall(FillLocalScalar(user, numerator, 0.0));
    PetscCall(FillLocalScalar(user, denominator, 0.0));
    PetscCall(FillLocalScalar(user, ratio, 0.0));
    PetscCall(SetInteriorScalar(user, numerator, 2.0));
    PetscCall(SetInteriorScalar(user, denominator, 1.0));
 
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvSpatialRatioAverage(user, &plan, numerator, denominator, NULL, all,
                                        PETSC_COMM_WORLD, ratio, NULL));
    PetscCall(ReadLocalScalar(user, ratio, 2, 2, 2, &clean));
 
    /* Poison both layout boundary planes in every direction. */
    PetscCall(DMDAVecGetArray(user->da, numerator, &num));
    for (PetscInt k = 0; k < info.mz; ++k)
    for (PetscInt j = 0; j < info.my; ++j)
    for (PetscInt i = 0; i < info.mx; ++i) {
        const PetscBool duplicate = (PetscBool)(i == 0 || j == 0 || k == 0 ||
                                                i == info.mx - 1 || j == info.my - 1 || k == info.mz - 1);
 
        if (duplicate) num[k][j][i] = 1.0e6;
    }
    PetscCall(DMDAVecRestoreArray(user->da, numerator, &num));
 
    PetscCall(SpatialTargetPlanCreate(user, FIELD_ID_CS,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvSpatialRatioAverage(user, &plan, numerator, denominator, NULL, all,
                                        PETSC_COMM_WORLD, ratio, NULL));
    PetscCall(ReadLocalScalar(user, ratio, 2, 2, 2, &poisoned));
 
    PetscCall(PicurvAssertRealNear(2.0, clean, 1.0e-13,
                                   "the interior average should be the interior value"));
    PetscCall(PicurvAssertRealNear(clean, poisoned, 1.0e-13,
                                   "the duplicate planes must not contribute to a spatial average"));
 
    PetscCall(VecDestroy(&numerator));
    PetscCall(VecDestroy(&denominator));
    PetscCall(VecDestroy(&ratio));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/** @brief Tests that homogeneous averaging falls back to the block's periodic axes. */
static PetscErrorCode TestHomogeneousAveragingDerivesPeriodicAxes(void)
{
    SimCtx   *simCtx = NULL;
    UserCtx  *user = NULL;
    PetscBool direction[3];
 
    PetscFunctionBeginUser;
    /* Periodic in xi and zeta only: the channel arrangement. The fixture sets the
       resolved periodicity flags, which are what the closure reads; the boundary faces
       are left alone here so the test cannot pass by consulting them instead. */
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 6, 6, 6,
                                                         PETSC_TRUE, PETSC_FALSE, PETSC_TRUE));
    simCtx->les_config.averaging_mode = LES_AVERAGING_HOMOGENEOUS;
 
    PetscCall(ResolveLESAveragingDirections(user, direction));
    PetscCall(PicurvAssertBool(direction[0], "a periodic xi axis should be averaged over"));
    PetscCall(PicurvAssertBool((PetscBool)(!direction[1]),
                               "a wall-normal axis should keep its own coefficient profile"));
    PetscCall(PicurvAssertBool(direction[2], "a periodic zeta axis should be averaged over"));
 
    /* An explicit selection overrides the periodic default. */
    simCtx->les_config.averaging_direction[0] = PETSC_TRUE;
    simCtx->les_config.averaging_direction[1] = PETSC_FALSE;
    simCtx->les_config.averaging_direction[2] = PETSC_FALSE;
    PetscCall(ResolveLESAveragingDirections(user, direction));
    PetscCall(PicurvAssertBool(direction[0], "an explicit direction should be honoured"));
    PetscCall(PicurvAssertBool((PetscBool)(!direction[2]),
                               "an explicit selection should replace the periodic default"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/** @brief Tests that local and global averaging ignore the configured direction list. */
static PetscErrorCode TestAveragingModesSelectTheirOwnDirections(void)
{
    SimCtx   *simCtx = NULL;
    UserCtx  *user = NULL;
    PetscBool direction[3];
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->les_config.averaging_direction[0] = PETSC_TRUE;
 
    simCtx->les_config.averaging_mode = LES_AVERAGING_LOCAL;
    PetscCall(ResolveLESAveragingDirections(user, direction));
    PetscCall(PicurvAssertBool((PetscBool)(!direction[0] && !direction[1] && !direction[2]),
                               "local averaging should select no direction at all"));
 
    simCtx->les_config.averaging_mode = LES_AVERAGING_GLOBAL;
    PetscCall(ResolveLESAveragingDirections(user, direction));
    PetscCall(PicurvAssertBool((PetscBool)(direction[0] && direction[1] && direction[2]),
                               "global averaging should select every direction"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the constant model builds its viscosity without a coefficient field.
 *
 * `UserCtx::CS` and `UserCtx::lCs` are deliberately left unallocated. If the constant
 * path ever reaches for them again this case crashes rather than silently reintroducing
 * a field of one repeated number.
 */
static PetscErrorCode TestConstantModelNeedsNoCoefficientField(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal ***nu_t = NULL;
    const PetscReal constant_cs = 0.2;
    /* u = x gives |S| = sqrt(2); with Aj = 1 the filter width is 1. */
    const PetscReal expected = constant_cs * constant_cs * PetscSqrtReal(2.0);
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
    simCtx->les = CONSTANT_SMAGORINSKY;
    simCtx->les_config.constant_cs = constant_cs;
 
    PetscCall(DMCreateGlobalVector(user->da, &user->Nu_t));
    PetscCall(DMCreateLocalVector(user->da, &user->lNu_t));
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(VecSet(user->Nu_t, 0.0));
 
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
    for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k)
    for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j)
    for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
        ucat[k][j][i].x = (PetscReal)i;
        ucat[k][j][i].y = 0.0;
        ucat[k][j][i].z = 0.0;
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Aj, INSERT_VALUES, user->lAj));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Aj, INSERT_VALUES, user->lAj));
 
    PetscCall(PicurvAssertBool((PetscBool)(user->lCs == NULL),
                               "the constant model should allocate no coefficient field"));
    PetscCall(ComputeEddyViscosityLES(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->Nu_t, &nu_t));
    PetscCall(PicurvAssertRealNear(expected, nu_t[2][2][2], 1.0e-9,
                                   "the constant model should read its coefficient from configuration"));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nu_t, &nu_t));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/** @brief Tests that the dynamic procedure refuses to run for the constant model. */
static PetscErrorCode TestDynamicProcedureRejectsConstantModel(void)
{
    SimCtx        *simCtx = NULL;
    UserCtx       *user = NULL;
    PetscErrorCode call_ierr;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
    simCtx->les = CONSTANT_SMAGORINSKY;
    simCtx->step = 5;
    simCtx->StartStep = 0;
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    call_ierr = ComputeSmagorinskyConstant(user);
    PetscCall(PetscPopErrorHandler());
 
    PetscCall(PicurvAssertBool((PetscBool)(call_ierr != 0),
                               "the dynamic procedure should reject a model that has no "
                               "coefficient field to fill"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the dynamic procedure returns a zero coefficient on uniform flow.
 *
 * Exercises the whole path: the strain precompute and its halo, the stencil gather,
 * every test filter, both contractions, the averaging, and the limiting. A uniform
 * field has no strain and no stress between the two filter widths, so the Germano
 * identity has nothing to fit and the coefficient must come out exactly zero.
 */
static PetscErrorCode TestDynamicProcedureVanishesOnUniformFlow(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal ***coefficient = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->les = DYNAMIC_SMAGORINSKY;
    simCtx->step = 5;
    simCtx->StartStep = 0;
 
    PetscCall(DMCreateGlobalVector(user->da, &user->CS));
    PetscCall(DMCreateLocalVector(user->da, &user->lCs));
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(VecSet(user->lAj, 1.0));
    PetscCall(VecSet(user->lNvert, 0.0));
 
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
    for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k)
    for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j)
    for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
        ucat[k][j][i].x = 1.25;
        ucat[k][j][i].y = -0.5;
        ucat[k][j][i].z = 3.0;
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
 
    PetscCall(ComputeSmagorinskyConstant(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->CS, &coefficient));
    for (PetscInt k = 1; k < 6; ++k)
    for (PetscInt j = 1; j < 6; ++j)
    for (PetscInt i = 1; i < 6; ++i) {
        PetscCall(PicurvAssertRealNear(0.0, coefficient[k][j][i], 1.0e-20,
                                       "a uniform field should give the dynamic procedure "
                                       "nothing to fit"));
    }
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->CS, &coefficient));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests the dynamic procedure end to end on a periodic sheared field.
 *
 * Runs the full procedure on a block that is periodic in all three directions, where
 * the strain precompute has to reach across the wrap correctly, and asks for global
 * averaging. Two things must hold: every cell receives the same finite coefficient,
 * because a global average is one number for the block; and the coefficient stays
 * inside the configured ceiling.
 */
static PetscErrorCode TestDynamicProcedureGlobalAverageIsUniform(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal ***coefficient = NULL;
    PetscReal    reference = 0.0;
    const FieldId cell_fields[] = {FIELD_ID_UCAT};
    const PetscReal wavenumber = 2.0 * PETSC_PI / 5.0; /* five distinct cells per direction */
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 6, 6, 6,
                                                         PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    DeclarePeriodicFaces(user, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE);
    simCtx->les = DYNAMIC_SMAGORINSKY;
    simCtx->step = 5;
    simCtx->StartStep = 0;
    simCtx->les_config.averaging_mode = LES_AVERAGING_GLOBAL;
 
    PetscCall(DMCreateGlobalVector(user->da, &user->CS));
    PetscCall(DMCreateLocalVector(user->da, &user->lCs));
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(VecSet(user->lAj, 1.0));
    PetscCall(VecSet(user->lNvert, 0.0));
 
    /* A divergence-free periodic shear: enough structure that both contractions are
       nonzero, and smooth enough that the coefficient stays in range. */
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
    for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k)
    for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j)
    for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
        ucat[k][j][i].x =  PetscSinReal(wavenumber * (PetscReal)j);
        ucat[k][j][i].y =  PetscSinReal(wavenumber * (PetscReal)k);
        ucat[k][j][i].z =  PetscSinReal(wavenumber * (PetscReal)i);
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
    PetscCall(SynchronizePeriodicCellFields(user, 1, cell_fields));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCAT));
 
    PetscCall(ComputeSmagorinskyConstant(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->CS, &coefficient));
    reference = coefficient[1][1][1];
    PetscCall(PicurvAssertBool((PetscBool)(PetscIsNormalReal(reference) || reference == 0.0),
                               "the dynamic coefficient must be finite"));
    PetscCall(PicurvAssertBool((PetscBool)(reference >= 0.0 &&
                                           reference <= simCtx->les_config.max_cs *
                                                        simCtx->les_config.max_cs),
                               "clamping should keep the coefficient inside its ceiling"));
    for (PetscInt k = 1; k < 6; ++k)
    for (PetscInt j = 1; j < 6; ++j)
    for (PetscInt i = 1; i < 6; ++i) {
        PetscCall(PicurvAssertRealNear(reference, coefficient[k][j][i], 1.0e-13,
                                       "global averaging should give the block one coefficient"));
    }
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->CS, &coefficient));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Runs the dynamic procedure with the Simpson kernel on a linear field and returns
 *        the raw coefficient at the block centre.
 * @param[in]  width_ratio Configured test-to-grid width ratio.
 * @param[out] coefficient The unclipped dynamic coefficient at cell (4, 4, 4).
 */
static PetscErrorCode SimpsonCoefficientOnLinearField(PetscReal width_ratio, PetscReal *coefficient)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal ***cs = NULL;
    const PetscReal gamma = 0.1;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 8, 8));
    simCtx->les = DYNAMIC_SMAGORINSKY;
    simCtx->step = 5;
    simCtx->StartStep = 0;
    simCtx->les_config.test_filter_kernel = LES_TEST_FILTER_SIMPSON_IK;
    simCtx->les_config.test_filter_width_ratio = width_ratio;
    simCtx->les_config.averaging_mode = LES_AVERAGING_LOCAL;
    simCtx->les_config.clip_mode = LES_CLIP_NONE;
 
    PetscCall(DMCreateGlobalVector(user->da, &user->CS));
    PetscCall(DMCreateLocalVector(user->da, &user->lCs));
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(VecSet(user->lAj, 1.0));
    PetscCall(VecSet(user->lNvert, 0.0));
 
    /* A divergence-free linear strain: the strain and both filtered strain terms are
       uniform, so the model tensor scales exactly with (alpha - 1), while the Leonard
       stress is the stencil's second moment along i and does not depend on alpha. */
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
    for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k)
    for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j)
    for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
        ucat[k][j][i].x = -gamma * (PetscReal)i;
        ucat[k][j][i].y =  gamma * (PetscReal)j;
        ucat[k][j][i].z =  0.0;
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
 
    PetscCall(ComputeSmagorinskyConstant(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->CS, &cs));
    *coefficient = cs[4][4][4];
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->CS, &cs));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the Simpson kernel's width ratio applies to the two directions it filters.
 *
 * Simpson filters i and k only, so the cube-root test-to-grid width ratio is
 * ratio^(2/3) and alpha is ratio^(4/3), not ratio^2. On a linear field the coefficient
 * is proportional to 1 / (alpha - 1), so two ratios fix the exponent without depending
 * on any other convention in the procedure.
 */
static PetscErrorCode TestDynamicProcedureSimpsonRatioCountsFilteredDirections(void)
{
    PetscReal coefficient_2 = 0.0, coefficient_4 = 0.0;
    const PetscReal alpha_2 = PetscPowReal(2.0, 4.0 / 3.0);
    const PetscReal alpha_4 = PetscPowReal(4.0, 4.0 / 3.0);
 
    PetscFunctionBeginUser;
    PetscCall(SimpsonCoefficientOnLinearField(2.0, &coefficient_2));
    PetscCall(SimpsonCoefficientOnLinearField(4.0, &coefficient_4));
    PetscCall(PicurvAssertBool((PetscBool)(coefficient_2 > 0.0 && coefficient_4 > 0.0),
                               "a linear strain should give a positive Simpson coefficient"));
    PetscCall(PicurvAssertRealNear((alpha_4 - 1.0) / (alpha_2 - 1.0), coefficient_2 / coefficient_4,
                                   1.0e-10,
                                   "the Simpson coefficient should scale as 1/(ratio^(4/3) - 1)"));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Entry point for the LES closure suite.
 */
int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"sym-tensor-algebra", TestSymTensorAlgebra},
        {"strain-rate-from-gradients", TestStrainRateFromGradients},
        {"filter-width-models-separate-on-stretched-cell", TestFilterWidthModelsSeparateOnStretchedCell},
        {"filter-width-is-independent-of-cell-orientation", TestFilterWidthIsIndependentOfCellOrientation},
        {"scotti-filter-width-matches-closed-form", TestScottiFilterWidthMatchesClosedForm},
        {"wale-and-vreman-kernels", TestWALEAndVremanKernels},
        {"leonard-stress-vanishes-on-uniform-flow", TestLeonardStressVanishesOnUniformFlow},
        {"germano-model-tensor-on-constant-strain", TestGermanoModelTensorOnConstantStrain},
        {"germano-model-tensor-uses-filtered-product", TestGermanoModelTensorUsesFilteredProduct},
        {"clip-model-coefficient-modes", TestClipModelCoefficientModes},
        {"eddy-viscosity-floor-bounds-total-viscosity", TestEddyViscosityFloorBoundsTotalViscosity},
        {"subgrid-kinetic-energy", TestSubgridKineticEnergy},
        {"average-ratio-local-is-pointwise", TestAverageRatioLocalIsPointwise},
        {"average-ratio-divides-summed-fields", TestAverageRatioDividesSummedFields},
        {"average-ratio-retains-unaveraged-direction", TestAverageRatioRetainsUnaveragedDirection},
        {"average-ratio-ignores-periodic-duplicate-planes", TestAverageRatioIgnoresPeriodicDuplicatePlanes},
        {"homogeneous-averaging-derives-periodic-axes", TestHomogeneousAveragingDerivesPeriodicAxes},
        {"averaging-modes-select-their-own-directions", TestAveragingModesSelectTheirOwnDirections},
        {"constant-model-needs-no-coefficient-field", TestConstantModelNeedsNoCoefficientField},
        {"dynamic-procedure-rejects-constant-model", TestDynamicProcedureRejectsConstantModel},
        {"dynamic-procedure-vanishes-on-uniform-flow", TestDynamicProcedureVanishesOnUniformFlow},
        {"dynamic-procedure-global-average-is-uniform", TestDynamicProcedureGlobalAverageIsUniform},
        {"dynamic-procedure-simpson-ratio-counts-filtered-directions", TestDynamicProcedureSimpsonRatioCountsFilteredDirections},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv LES closure tests");
    if (ierr) {
        return (int)ierr;
    }
 
    ierr = PicurvRunTests("unit-les", cases, sizeof(cases) / sizeof(cases[0]));
    if (ierr) {
        PetscFinalize();
        return (int)ierr;
    }
 
    ierr = PetscFinalize();
    return (int)ierr;
}