test_solver_kernels.c File Reference

C unit tests for solver-side analytical and LES helper kernels. More...

#include "test_support.h"
#include "AnalyticalSolutions.h"
#include "BodyForces.h"
#include "Filter.h"
#include "les.h"
#include "solvers.h"
#include "momentumsolvers.h"

Include dependency graph for test_solver_kernels.c:

Go to the source code of this file.

Functions
static PetscErrorCode	TestLESTestFilterPaths (void)
	Tests LES test-filter helper paths for representative cases.
static PetscErrorCode	TestAnalyticalGeometrySelection (void)
	Tests analytical geometry selection for supported analytical solutions.
static PetscErrorCode	TestAnalyticalScalarVerificationHelpers (void)
	Tests analytical scalar verification helper routines.
static PetscErrorCode	TestAnalyticalSolutionEngineDispatch (void)
	Tests analytical solution engine ZERO_FLOW, UNIFORM_FLOW, and unknown-type dispatch.
static PetscErrorCode	TestAnalyticalSolutionEngineTaylorGreenSamples (void)
	Tests exact Taylor-Green samples on selected Eulerian interior and boundary points.
static PetscErrorCode	TestAnalyticalSolutionForParticlesDispatch (void)
	Tests particle analytical-solution dispatch for TGV3D, UNIFORM_FLOW, and non-analytical no-op paths.
static PetscErrorCode	TestComputeEddyViscosityLESDeterministicField (void)
	Tests deterministic LES eddy-viscosity computation on a linear velocity field.
static PetscErrorCode	TestFlowSolverRejectsUnsupportedMomentumSolverType (void)
	Tests FlowSolver guardrails for unsupported momentum solver selections.
static PetscErrorCode	TestDrivenChannelFlowSource (void)
	Tests driven-channel flow source-term evaluation.
static void	MomSetAllBC (UserCtx *user, BCType t)
	Sets the mathematical BC type on all six faces of a test UserCtx.
static PetscErrorCode	MomFillLocalCmpnts (DM fda, Vec lvec, PetscReal x, PetscReal y, PetscReal z)
	Fills a local Cmpnts vector (ghosts included) with a uniform vector value.
static PetscErrorCode	MomFillLocalScalar (DM da, Vec lvec, PetscReal v)
	Fills a local scalar vector (ghosts included) with a uniform value.
static PetscErrorCode	MomSetLocalScalarCell (DM da, Vec lvec, PetscInt ci, PetscInt cj, PetscInt ck, PetscReal v)
	Sets a single local-scalar cell value (e.g.
static PetscErrorCode	TestMomentumBDFCoefficient (void)
static PetscErrorCode	MomMakeUnitGridP (SimCtx **simCtx, UserCtx **user, PetscInt n, PetscBool px, PetscBool py, PetscBool pz, BCType bc)
static PetscErrorCode	MomMakeUnitGrid (SimCtx **simCtx, UserCtx **user, PetscInt n, BCType bc)
static PetscErrorCode	TestMomentumStabilityCenteredConvection (void)
static PetscErrorCode	TestMomentumStabilityQuickInterior (void)
static PetscErrorCode	TestMomentumStabilityQuickBoundary (void)
static PetscErrorCode	TestMomentumStabilityViscousLongitudinal (void)
static PetscErrorCode	TestMomentumStabilityViscousScalesWithNu (void)
static PetscErrorCode	TestMomentumStabilityLESEddyViscosity (void)
static PetscErrorCode	TestMomentumStabilityOneSidedViscous (void)
static PetscErrorCode	TestMomentumStabilityReadOnly (void)
static PetscErrorCode	MomSetCartesianMetrics (UserCtx *user, PetscReal dx, PetscReal dy, PetscReal dz)
	Sets uniform Cartesian metrics for spacing (dx,dy,dz) on all cell and face arrays.
static PetscErrorCode	MomFillUcontXRamp (DM fda, Vec lvec, PetscReal slope)
	Fills lUcont.x as a linear ramp slope*i (nonzero discrete contravariant divergence).
static PetscErrorCode	MomFillUcatShearY (DM fda, Vec lvec, PetscReal gamma)
	Fills lUcat as a linear shear u=(gamma*j,0,0) so |grad u|_inf = gamma.
static PetscErrorCode	TestMomentumStabilityAnisotropicViscous (void)
static PetscErrorCode	TestMomentumStabilityDirectionalQuick (void)
static PetscErrorCode	TestMomentumStabilityOneSidedDirections (void)
static PetscErrorCode	TestMomentumStabilityActiveRowMask (void)
static PetscErrorCode	TestMomentumStabilityValidation (void)
static PetscErrorCode	TestMomentumStabilityCandidateBCDivergence (void)
static PetscErrorCode	TestMomentumStabilityCandidateDShear (void)
static PetscErrorCode	TestMomentumStabilityWallSuppression (void)
static PetscErrorCode	TestMomentumStabilityReadOnlyExact (void)
static PetscErrorCode	TestMomentumActiveRowsHelper (void)
static PetscErrorCode	TestMomentumStabilityErrorReentry (void)
static PetscErrorCode	TestMomentumBDFResidualCombination (void)
int	main (int argc, char **argv)
	Runs the unit-solver PETSc test binary.

Detailed Description

C unit tests for solver-side analytical and LES helper kernels.

Definition in file test_solver_kernels.c.

Function Documentation

TestLESTestFilterPaths()

static PetscErrorCode TestLESTestFilterPaths ( void  )

static

Tests LES test-filter helper paths for representative cases.

Definition at line 18 of file test_solver_kernels.c.

{
    SimCtx simCtx;
    double values[3][3][3];
    double weights[3][3][3];
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&simCtx, sizeof(simCtx)));
    for (PetscInt k = 0; k < 3; ++k) {
        for (PetscInt j = 0; j < 3; ++j) {
            for (PetscInt i = 0; i < 3; ++i) {
                values[k][j][i] = 2.0;
                weights[k][j][i] = 1.0;
            }
        }
    }
 
    simCtx.testfilter_ik = 1;
    PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
                                   "Simpson-rule filter should preserve a constant field"));
 
    simCtx.testfilter_ik = 0;
    PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
                                   "box filter should preserve a constant field"));
 
    for (PetscInt k = 0; k < 3; ++k) {
        for (PetscInt j = 0; j < 3; ++j) {
            for (PetscInt i = 0; i < 3; ++i) {
                weights[k][j][i] = 0.0;
            }
        }
    }
    PetscCall(PicurvAssertRealNear(0.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
                                   "box filter should return zero when all weights are zero"));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAnalyticalGeometrySelection()

static PetscErrorCode TestAnalyticalGeometrySelection ( void  )

static

Tests analytical geometry selection for supported analytical solutions.

Definition at line 58 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscErrorCode ierr_grid = 0;
    PetscErrorCode ierr_non_square = 0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 8, 8));
    PetscCall(PicurvAssertBool(AnalyticalTypeRequiresCustomGeometry("TGV3D"),
                               "TGV3D should require custom geometry"));
    PetscCall(PicurvAssertBool((PetscBool)!AnalyticalTypeRequiresCustomGeometry("ZERO_FLOW"),
                               "ZERO_FLOW should not require custom geometry"));
 
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    ierr_grid = SetAnalyticalGridInfo(user);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(ierr_grid != 0),
                               "SetAnalyticalGridInfo should reject analytical types without custom geometry"));
 
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "TGV3D", sizeof(simCtx->AnalyticalSolutionType)));
    simCtx->block_number = 2;
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    ierr_non_square = SetAnalyticalGridInfo(user);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(ierr_non_square != 0),
                               "TGV3D multi-block setup should reject non-square block counts"));
 
    simCtx->block_number = 1;
    PetscCall(SetAnalyticalGridInfo(user));
    PetscCall(PicurvAssertRealNear(0.0, user->Min_X, 1.0e-12, "TGV3D single-block xmin"));
    PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_X, 1.0e-12, "TGV3D single-block xmax"));
 
    simCtx->block_number = 4;
    user->_this = 3;
    PetscCall(SetAnalyticalGridInfo(user));
    PetscCall(PicurvAssertRealNear(PETSC_PI, user->Min_X, 1.0e-12, "TGV3D multi-block xmin should reflect the block column"));
    PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_X, 1.0e-12, "TGV3D multi-block xmax should reflect the block column"));
    PetscCall(PicurvAssertRealNear(PETSC_PI, user->Min_Y, 1.0e-12, "TGV3D multi-block ymin should reflect the block row"));
    PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_Y, 1.0e-12, "TGV3D multi-block ymax should reflect the block row"));
    PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_Z, 1.0e-12, "TGV3D multi-block zmax should span the full domain"));
 
    simCtx->block_number = 1;
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAnalyticalScalarVerificationHelpers()

static PetscErrorCode TestAnalyticalScalarVerificationHelpers ( void  )

static

Tests analytical scalar verification helper routines.

Definition at line 109 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec target = NULL;
    PetscReal value = 0.0;
    PetscReal ***target_arr = NULL;
    PetscReal *positions = NULL;
    PetscReal *psi = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
    simCtx->verificationScalar.enabled = PETSC_TRUE;
    PetscCall(PetscStrncpy(simCtx->verificationScalar.mode,
                           "analytical",
                           sizeof(simCtx->verificationScalar.mode)));
 
    PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
                           "CONSTANT",
                           sizeof(simCtx->verificationScalar.profile)));
    simCtx->verificationScalar.value = 2.5;
    PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.2, 0.3, 0.4, 0.0, &value));
    PetscCall(PicurvAssertRealNear(2.5, value, 1.0e-12,
                                   "constant scalar profile should evaluate to the configured value"));
 
    PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
    positions[0] = 0.1; positions[1] = 0.2; positions[2] = 0.3;
    positions[3] = 0.8; positions[4] = 0.6; positions[5] = 0.4;
    PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
 
    PetscCall(SetAnalyticalScalarFieldOnParticles(user, "Psi"));
    PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    PetscCall(PicurvAssertRealNear(2.5, psi[0], 1.0e-12,
                                   "SetAnalyticalScalarFieldOnParticles should overwrite the first particle scalar"));
    PetscCall(PicurvAssertRealNear(2.5, psi[1], 1.0e-12,
                                   "SetAnalyticalScalarFieldOnParticles should overwrite the second particle scalar"));
    PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
 
    PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
                           "LINEAR_X",
                           sizeof(simCtx->verificationScalar.profile)));
    simCtx->verificationScalar.phi0 = 1.0;
    simCtx->verificationScalar.slope_x = 2.0;
    PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.25, 0.0, 0.0, 0.0, &value));
    PetscCall(PicurvAssertRealNear(1.5, value, 1.0e-12,
                                   "linear-x scalar profile should evaluate phi0 + slope_x * x"));
 
    PetscCall(VecDuplicate(user->Psi, &target));
    PetscCall(SetAnalyticalScalarFieldAtCellCenters(user, target));
    PetscCall(DMDAVecGetArrayRead(user->da, target, &target_arr));
    PetscCall(PicurvAssertRealNear(1.25, target_arr[1][1][1], 1.0e-12,
                                   "cell-center scalar fill should use the physical x center coordinate"));
    PetscCall(DMDAVecRestoreArrayRead(user->da, target, &target_arr));
    PetscCall(VecDestroy(&target));
 
    PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
                           "SIN_PRODUCT",
                           sizeof(simCtx->verificationScalar.profile)));
    simCtx->verificationScalar.amplitude = 3.0;
    simCtx->verificationScalar.kx = PETSC_PI;
    simCtx->verificationScalar.ky = PETSC_PI;
    simCtx->verificationScalar.kz = PETSC_PI;
    PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.5, 0.5, 0.5, 0.0, &value));
    PetscCall(PicurvAssertRealNear(3.0, value, 1.0e-12,
                                   "sin-product scalar profile should peak at pi/2 in each coordinate"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAnalyticalSolutionEngineDispatch()

static PetscErrorCode TestAnalyticalSolutionEngineDispatch ( void  )

static

Tests analytical solution engine ZERO_FLOW, UNIFORM_FLOW, and unknown-type dispatch.

Definition at line 183 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscErrorCode ierr_unknown = 0;
    Cmpnts ***ucat = NULL;
    Cmpnts ***ubcs = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(VecSet(user->Ucat, 3.0));
    PetscCall(VecSet(user->P, 5.0));
    PetscCall(VecSet(user->Bcs.Ubcs, 7.0));
 
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
    PetscCall(AnalyticalSolutionEngine(simCtx));
    PetscCall(PicurvAssertVecConstant(user->Ucat, 0.0, 1.0e-12, "ZERO_FLOW should zero the Eulerian velocity field"));
    PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "ZERO_FLOW should zero the pressure field"));
    PetscCall(PicurvAssertVecConstant(user->Bcs.Ubcs, 0.0, 1.0e-12, "ZERO_FLOW should zero boundary-condition velocity data"));
 
    /* UNIFORM_FLOW now works in curvilinear form (sets Ucont via metric dot products,
       derives Ucat via Contra2Cart).  The test grid needs identity metrics so the
       transformation is invertible and Ucat recovers the physical velocity exactly. */
    {
        Cmpnts ***l_csi, ***l_eta, ***l_zet;
        DMDALocalInfo linfo;
        PetscCall(DMDAGetLocalInfo(user->fda, &linfo));
        PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
        PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
        PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
        for (PetscInt k = linfo.zs; k < linfo.zs + linfo.zm; k++)
            for (PetscInt j = linfo.ys; j < linfo.ys + linfo.ym; j++)
                for (PetscInt i = linfo.xs; i < linfo.xs + linfo.xm; i++) {
                    l_csi[k][j][i] = (Cmpnts){1.0, 0.0, 0.0};
                    l_eta[k][j][i] = (Cmpnts){0.0, 1.0, 0.0};
                    l_zet[k][j][i] = (Cmpnts){0.0, 0.0, 1.0};
                }
        PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
        PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
        PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
    }
 
    simCtx->AnalyticalUniformVelocity.x = 1.25;
    simCtx->AnalyticalUniformVelocity.y = -0.5;
    simCtx->AnalyticalUniformVelocity.z = 0.75;
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
    PetscCall(AnalyticalSolutionEngine(simCtx));
    PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "UNIFORM_FLOW should keep the pressure field zero"));
 
    Cmpnts ***ucont = NULL;
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(PicurvAssertRealNear(1.25, ucat[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should impose the configured x velocity"));
    PetscCall(PicurvAssertRealNear(-0.5, ucat[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should impose the configured y velocity"));
    PetscCall(PicurvAssertRealNear(0.75, ucat[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should impose the configured z velocity"));
    PetscCall(PicurvAssertRealNear(1.25, ucont[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should set contravariant x flux (identity metric)"));
    PetscCall(PicurvAssertRealNear(-0.5, ucont[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should set contravariant y flux (identity metric)"));
    PetscCall(PicurvAssertRealNear(0.75, ucont[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should set contravariant z flux (identity metric)"));
    PetscCall(PicurvAssertRealNear(1.25, ubcs[0][1][1].x, 1.0e-12, "UNIFORM_FLOW should populate boundary x velocity"));
    PetscCall(PicurvAssertRealNear(-0.5, ubcs[0][1][1].y, 1.0e-12, "UNIFORM_FLOW should populate boundary y velocity"));
    PetscCall(PicurvAssertRealNear(0.75, ubcs[0][1][1].z, 1.0e-12, "UNIFORM_FLOW should populate boundary z velocity"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
 
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "NOT_A_REAL_ANALYTICAL_TYPE", sizeof(simCtx->AnalyticalSolutionType)));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    ierr_unknown = AnalyticalSolutionEngine(simCtx);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(ierr_unknown != 0),
                               "AnalyticalSolutionEngine should reject unknown analytical type strings"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAnalyticalSolutionEngineTaylorGreenSamples()

static PetscErrorCode TestAnalyticalSolutionEngineTaylorGreenSamples ( void  )

static

Tests exact Taylor-Green samples on selected Eulerian interior and boundary points.

Definition at line 262 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Cmpnts ***cent = NULL;
    Cmpnts ***cent_x = NULL;
    Cmpnts ***cent_y = NULL;
    Cmpnts ***cent_z = NULL;
    Cmpnts ***ucat = NULL;
    Cmpnts ***ubcs = NULL;
    PetscReal ***p = NULL;
    const PetscReal vel_decay = PetscExpReal(-0.5);
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    simCtx->ren = 2.0;
    simCtx->ti = 0.5;
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "TGV3D", sizeof(simCtx->AnalyticalSolutionType)));
 
    PetscCall(DMDAVecGetArray(user->fda, user->Cent, &cent));
    PetscCall(DMDAVecGetArray(user->fda, user->lCentx, &cent_x));
    PetscCall(DMDAVecGetArray(user->fda, user->lCenty, &cent_y));
    PetscCall(DMDAVecGetArray(user->fda, user->lCentz, &cent_z));
    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) {
                cent[k][j][i].x = 0.0;
                cent[k][j][i].y = 0.0;
                cent[k][j][i].z = 0.0;
                cent_x[k][j][i] = cent[k][j][i];
                cent_y[k][j][i] = cent[k][j][i];
                cent_z[k][j][i] = cent[k][j][i];
            }
        }
    }
    cent[1][1][1].x = 0.5 * PETSC_PI;
    cent[1][1][1].y = 0.0;
    cent[1][1][1].z = 0.0;
    cent[1][2][1].x = 0.0;
    cent[1][2][1].y = 0.5 * PETSC_PI;
    cent[1][2][1].z = 0.0;
    cent_z[0][1][1].x = 0.5 * PETSC_PI;
    cent_z[0][1][1].y = 0.0;
    cent_z[0][1][1].z = 0.0;
    PetscCall(DMDAVecRestoreArray(user->fda, user->lCentz, &cent_z));
    PetscCall(DMDAVecRestoreArray(user->fda, user->lCenty, &cent_y));
    PetscCall(DMDAVecRestoreArray(user->fda, user->lCentx, &cent_x));
    PetscCall(DMDAVecRestoreArray(user->fda, user->Cent, &cent));
 
    PetscCall(AnalyticalSolutionEngine(simCtx));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(DMDAVecGetArrayRead(user->da, user->P, &p));
    PetscCall(PicurvAssertRealNear(vel_decay, ucat[1][1][1].x, 1.0e-12, "TGV sample should set the expected interior x velocity"));
    PetscCall(PicurvAssertRealNear(0.0, ucat[1][1][1].y, 1.0e-12, "TGV sample should keep the paired interior y velocity at zero"));
    PetscCall(PicurvAssertRealNear(-vel_decay, ucat[1][2][1].y, 1.0e-12, "TGV sample should set the expected interior y velocity"));
    PetscCall(PicurvAssertRealNear(0.0, p[1][1][1], 1.0e-12, "Chosen TGV sample should produce zero pressure"));
    PetscCall(PicurvAssertRealNear(vel_decay, ubcs[0][1][1].x, 1.0e-12, "TGV sample should set the expected boundary x velocity"));
    PetscCall(PicurvAssertRealNear(0.0, ubcs[0][1][1].y, 1.0e-12, "TGV boundary sample should keep the paired y velocity at zero"));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->P, &p));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAnalyticalSolutionForParticlesDispatch()

static PetscErrorCode TestAnalyticalSolutionForParticlesDispatch ( void  )

static

Tests particle analytical-solution dispatch for TGV3D, UNIFORM_FLOW, and non-analytical no-op paths.

Definition at line 333 of file test_solver_kernels.c.

{
    SimCtx simCtx;
    Vec tempVec = NULL;
    PetscReal *data = NULL;
    const PetscReal vel_decay = PetscExpReal(-0.5);
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&simCtx, sizeof(simCtx)));
    simCtx.ren = 2.0;
    simCtx.ti = 0.5;
    PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "TGV3D", sizeof(simCtx.AnalyticalSolutionType)));
 
    PetscCall(VecCreateSeq(PETSC_COMM_SELF, 6, &tempVec));
    PetscCall(VecGetArray(tempVec, &data));
    data[0] = 0.5 * PETSC_PI; data[1] = 0.0;          data[2] = 0.0;
    data[3] = 0.0;           data[4] = 0.5 * PETSC_PI; data[5] = 0.0;
    PetscCall(VecRestoreArray(tempVec, &data));
 
    PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
    PetscCall(VecGetArray(tempVec, &data));
    PetscCall(PicurvAssertRealNear(vel_decay, data[0], 1.0e-12,
                                   "TGV3D particle dispatch should populate the x velocity at x=pi/2"));
    PetscCall(PicurvAssertRealNear(0.0, data[1], 1.0e-12,
                                   "TGV3D particle dispatch should leave the first particle y velocity at zero"));
    PetscCall(PicurvAssertRealNear(0.0, data[2], 1.0e-12,
                                   "TGV3D particle dispatch should leave the first particle z velocity at zero"));
    PetscCall(PicurvAssertRealNear(0.0, data[3], 1.0e-12,
                                   "TGV3D particle dispatch should leave the second particle x velocity at zero"));
    PetscCall(PicurvAssertRealNear(-vel_decay, data[4], 1.0e-12,
                                   "TGV3D particle dispatch should populate the y velocity at y=pi/2"));
    PetscCall(PicurvAssertRealNear(0.0, data[5], 1.0e-12,
                                   "TGV3D particle dispatch should leave the second particle z velocity at zero"));
    PetscCall(VecRestoreArray(tempVec, &data));
 
    PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
    PetscCall(VecGetArray(tempVec, &data));
    data[0] = 3.0;
    data[1] = 4.0;
    data[2] = 5.0;
    PetscCall(VecRestoreArray(tempVec, &data));
    PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
    PetscCall(VecGetArray(tempVec, &data));
    PetscCall(PicurvAssertRealNear(3.0, data[0], 1.0e-12,
                                   "Non-TGV particle dispatch should leave the vector untouched"));
    PetscCall(PicurvAssertRealNear(4.0, data[1], 1.0e-12,
                                   "Non-TGV particle dispatch should preserve the y component"));
    PetscCall(PicurvAssertRealNear(5.0, data[2], 1.0e-12,
                                   "Non-TGV particle dispatch should preserve the z component"));
    PetscCall(VecRestoreArray(tempVec, &data));
 
    simCtx.AnalyticalUniformVelocity.x = 0.125;
    simCtx.AnalyticalUniformVelocity.y = -0.25;
    simCtx.AnalyticalUniformVelocity.z = 0.375;
    PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
    PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
    PetscCall(VecGetArray(tempVec, &data));
    PetscCall(PicurvAssertRealNear(0.125, data[0], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should populate the x velocity"));
    PetscCall(PicurvAssertRealNear(-0.25, data[1], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should populate the y velocity"));
    PetscCall(PicurvAssertRealNear(0.375, data[2], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should populate the z velocity"));
    PetscCall(PicurvAssertRealNear(0.125, data[3], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should use the same x velocity for each particle"));
    PetscCall(PicurvAssertRealNear(-0.25, data[4], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should use the same y velocity for each particle"));
    PetscCall(PicurvAssertRealNear(0.375, data[5], 1.0e-12,
                                   "UNIFORM_FLOW particle dispatch should use the same z velocity for each particle"));
    PetscCall(VecRestoreArray(tempVec, &data));
 
    PetscCall(VecDestroy(&tempVec));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestComputeEddyViscosityLESDeterministicField()

static PetscErrorCode TestComputeEddyViscosityLESDeterministicField ( void  )

static

Tests deterministic LES eddy-viscosity computation on a linear velocity field.

Definition at line 411 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Cmpnts ***ucat = NULL;
    PetscReal ***nu_t = NULL;
    const PetscReal expected_nu_t = 0.25 * PetscSqrtReal(2.0);
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
    PetscCall(DMCreateGlobalVector(user->da, &user->Nu_t));
    PetscCall(DMCreateLocalVector(user->da, &user->lNu_t));
    PetscCall(DMCreateGlobalVector(user->da, &user->CS));
    PetscCall(DMCreateLocalVector(user->da, &user->lCs));
 
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(VecSet(user->Nu_t, 0.0));
    PetscCall(VecSet(user->CS, 0.5));
    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(DMGlobalToLocalBegin(user->da, user->CS, INSERT_VALUES, user->lCs));
    PetscCall(DMGlobalToLocalEnd(user->da, user->CS, INSERT_VALUES, user->lCs));
 
    PetscCall(ComputeEddyViscosityLES(user));
    PetscCall(DMDAVecGetArrayRead(user->da, user->Nu_t, &nu_t));
    PetscCall(PicurvAssertRealNear(expected_nu_t, nu_t[2][2][2], 1.0e-6,
                                   "linear velocity field should yield deterministic LES eddy viscosity"));
    PetscCall(PicurvAssertRealNear(0.0, nu_t[0][2][2], 1.0e-12,
                                   "boundary cells should remain untouched by the interior LES loop"));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nu_t, &nu_t));
 
    PetscCall(VecDestroy(&user->CS));
    PetscCall(VecDestroy(&user->lCs));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestFlowSolverRejectsUnsupportedMomentumSolverType()

static PetscErrorCode TestFlowSolverRejectsUnsupportedMomentumSolverType ( void  )

static

Tests FlowSolver guardrails for unsupported momentum solver selections.

Definition at line 464 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscErrorCode ierr_flow = 0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    simCtx->mom_solver_type = (MomentumSolverType)999;
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    ierr_flow = FlowSolver(simCtx);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(ierr_flow != 0),
                               "FlowSolver should reject unsupported momentum solver selectors"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestDrivenChannelFlowSource()

static PetscErrorCode TestDrivenChannelFlowSource ( void  )

static

Tests driven-channel flow source-term evaluation.

Definition at line 487 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec rct = NULL;
    Cmpnts ***rct_arr = NULL;
    Cmpnts ***l_csi = NULL;
    Cmpnts ***l_eta = NULL;
    Cmpnts ***l_zet = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(VecDuplicate(user->Ucont, &rct));
    PetscCall(VecZeroEntries(rct));
 
    user->boundary_faces[BC_FACE_NEG_X].face_id = BC_FACE_NEG_X;
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = PERIODIC;
    user->boundary_faces[BC_FACE_NEG_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
    user->boundary_faces[BC_FACE_POS_X].face_id = BC_FACE_POS_X;
    user->boundary_faces[BC_FACE_POS_X].mathematical_type = PERIODIC;
    user->boundary_faces[BC_FACE_POS_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
    simCtx->bulkVelocityCorrection = 2.0;
    simCtx->dt = 1.0;
    simCtx->forceScalingFactor = 1.0;
    simCtx->drivingForceMagnitude = 0.0;
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
    PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
    PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
    l_csi[1][1][1].x = 1.0; l_csi[1][1][1].y = 0.0; l_csi[1][1][1].z = 0.0;
    l_eta[1][1][1].x = 0.0; l_eta[1][1][1].y = 1.0; l_eta[1][1][1].z = 0.0;
    l_zet[1][1][1].x = 0.0; l_zet[1][1][1].y = 0.0; l_zet[1][1][1].z = 1.0;
    PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
    PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
    PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
 
    PetscCall(ComputeDrivenChannelFlowSource(user, rct));
    PetscCall(DMDAVecGetArrayRead(user->fda, rct, &rct_arr));
    PetscCall(PicurvAssertRealNear(1.5, simCtx->drivingForceMagnitude, 1.0e-12,
                                   "driven flow source should update the controller magnitude"));
    PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].y, 1.0e-12,
                                   "driven flow source should leave the y component unchanged"));
    PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].z, 1.0e-12,
                                   "driven flow source should leave the z component unchanged"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, rct, &rct_arr));
 
    PetscCall(VecDestroy(&rct));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomSetAllBC()

static void MomSetAllBC ( UserCtx *  user, BCType  t  )

static

Sets the mathematical BC type on all six faces of a test UserCtx.

Parameters

user	Test user context.
t	BCType to assign to every face (e.g. PERIODIC, INLET).

Definition at line 551 of file test_solver_kernels.c.

{
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = t;
    user->boundary_faces[BC_FACE_POS_X].mathematical_type = t;
    user->boundary_faces[BC_FACE_NEG_Y].mathematical_type = t;
    user->boundary_faces[BC_FACE_POS_Y].mathematical_type = t;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = t;
    user->boundary_faces[BC_FACE_POS_Z].mathematical_type = t;
}

MomFillLocalCmpnts()

static PetscErrorCode MomFillLocalCmpnts ( DM  fda, Vec  lvec, PetscReal  x, PetscReal  y, PetscReal  z  )

static

Fills a local Cmpnts vector (ghosts included) with a uniform vector value.

Parameters

fda	Vector DM owning lvec.
lvec	Local Cmpnts vector to fill.
x,y,z	Uniform component values.

Returns

PetscErrorCode 0 on success.

Definition at line 568 of file test_solver_kernels.c.

{
    Cmpnts ***a;
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(fda, &gxs, &gys, &gzs, &gxm, &gym, &gzm));
    PetscCall(DMDAVecGetArray(fda, lvec, &a));
    for (PetscInt k = gzs; k < gzs+gzm; ++k)
        for (PetscInt j = gys; j < gys+gym; ++j)
            for (PetscInt i = gxs; i < gxs+gxm; ++i) { a[k][j][i].x = x; a[k][j][i].y = y; a[k][j][i].z = z; }
    PetscCall(DMDAVecRestoreArray(fda, lvec, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MomFillLocalScalar()

static PetscErrorCode MomFillLocalScalar ( DM  da, Vec  lvec, PetscReal  v  )

static

Fills a local scalar vector (ghosts included) with a uniform value.

Parameters

da	Scalar DM owning lvec.
lvec	Local scalar vector to fill.
v	Uniform value.

Returns

PetscErrorCode 0 on success.

Definition at line 589 of file test_solver_kernels.c.

{
    PetscReal ***a;
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gxm, &gym, &gzm));
    PetscCall(DMDAVecGetArray(da, lvec, &a));
    for (PetscInt k = gzs; k < gzs+gzm; ++k)
        for (PetscInt j = gys; j < gys+gym; ++j)
            for (PetscInt i = gxs; i < gxs+gxm; ++i) a[k][j][i] = v;
    PetscCall(DMDAVecRestoreArray(da, lvec, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MomSetLocalScalarCell()

static PetscErrorCode MomSetLocalScalarCell ( DM  da, Vec  lvec, PetscInt  ci, PetscInt  cj, PetscInt  ck, PetscReal  v  )

static

Sets a single local-scalar cell value (e.g.

to mark one nvert solid).

Parameters

da	Scalar DM owning lvec.
lvec	Local scalar vector.
ci,cj,ck	Cell indices to set.
v	Value to assign.

Returns

PetscErrorCode 0 on success.

Definition at line 611 of file test_solver_kernels.c.

{
    PetscReal ***a;
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gxm, &gym, &gzm));
    PetscCall(DMDAVecGetArray(da, lvec, &a));
    /* MPI-safe: only write if (ci,cj,ck) is within this rank's ghost-inclusive range. */
    if (ci >= gxs && ci < gxs+gxm && cj >= gys && cj < gys+gym && ck >= gzs && ck < gzs+gzm)
        a[ck][cj][ci] = v;
    PetscCall(DMDAVecRestoreArray(da, lvec, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestMomentumBDFCoefficient()

static PetscErrorCode TestMomentumBDFCoefficient ( void  )

static

Definition at line 626 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL;
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->StartStep = 0; simCtx->step = 1;
    PetscCall(PicurvAssertRealNear(1.0, MomentumBDFCoefficient(simCtx), 1e-12, "a0 BDF1 on cold-start step 1"));
    PetscCall(PicurvAssertBool((PetscBool)(!MomentumUsesBDF2(simCtx)), "step 1 is BDF1"));
    simCtx->step = 5;
    PetscCall(PicurvAssertRealNear(1.5, MomentumBDFCoefficient(simCtx), 1e-12, "a0 BDF2 on interior step"));
    PetscCall(PicurvAssertBool(MomentumUsesBDF2(simCtx), "step 5 (StartStep 0) is BDF2"));
    simCtx->StartStep = 5; simCtx->step = 5;
    PetscCall(PicurvAssertRealNear(1.0, MomentumBDFCoefficient(simCtx), 1e-12, "a0 BDF1 on restart step (ti==tistart)"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomMakeUnitGridP()

static PetscErrorCode MomMakeUnitGridP ( SimCtx **  simCtx, UserCtx **  user, PetscInt  n, PetscBool  px, PetscBool  py, PetscBool  pz, BCType  bc  )

static

Definition at line 646 of file test_solver_kernels.c.

{
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(simCtx, user, n, n, n, px, py, pz));
    (*user)->boundary_faces[BC_FACE_NEG_X].mathematical_type = px ? PERIODIC : bc;
    (*user)->boundary_faces[BC_FACE_POS_X].mathematical_type = px ? PERIODIC : bc;
    (*user)->boundary_faces[BC_FACE_NEG_Y].mathematical_type = py ? PERIODIC : bc;
    (*user)->boundary_faces[BC_FACE_POS_Y].mathematical_type = py ? PERIODIC : bc;
    (*user)->boundary_faces[BC_FACE_NEG_Z].mathematical_type = pz ? PERIODIC : bc;
    (*user)->boundary_faces[BC_FACE_POS_Z].mathematical_type = pz ? PERIODIC : bc;
    (*simCtx)->dt = 0.1; (*simCtx)->step = 1; (*simCtx)->StartStep = 0;   /* a0=1 -> lambda_t=10 */
    (*simCtx)->ren = 1.0; (*simCtx)->les = 0; (*simCtx)->rans = 0;
    (*simCtx)->central = 0; (*simCtx)->invicid = 0; (*simCtx)->block_number = 1;
    (*simCtx)->TwoD = 0; (*simCtx)->clark = 0;
    /* The minimal fixture does not allocate lNu_t (LES off by default); create a
       zeroed one so the estimator can read it when a test enables LES/RANS. */
    if (!(*user)->lNu_t) PetscCall(DMCreateLocalVector((*user)->da, &(*user)->lNu_t));
    PetscCall(MomFillLocalCmpnts((*user)->fda, (*user)->lUcont, 0.0, 0.0, 0.0));
    PetscCall(MomFillLocalCmpnts((*user)->fda, (*user)->lUcat,  0.0, 0.0, 0.0));
    PetscCall(MomFillLocalScalar((*user)->da, (*user)->lNvert, 0.0));
    PetscCall(MomFillLocalScalar((*user)->da, (*user)->lNu_t,  0.0));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomMakeUnitGrid()

static PetscErrorCode MomMakeUnitGrid ( SimCtx **  simCtx, UserCtx **  user, PetscInt  n, BCType  bc  )

static

Definition at line 672 of file test_solver_kernels.c.

{
    const PetscBool p = (PetscBool)(bc == PERIODIC);
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGridP(simCtx, user, n, p, p, p, bc));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityCenteredConvection()

static PetscErrorCode TestMomentumStabilityCenteredConvection ( void  )

static

Definition at line 681 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->central = 1; simCtx->invicid = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 2.0, 0.0, 0.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(10.0, rep.lambda_t, 1e-9, "lambda_t = a0/dt"));
    PetscCall(PicurvAssertRealNear(2.0,  rep.lambda_c, 1e-9, "centered conv f_c=1 (Aj*U^xi scaling)"));
    PetscCall(PicurvAssertRealNear(0.0,  rep.lambda_v, 1e-9, "inviscid: no viscous term"));
    PetscCall(PicurvAssertRealNear(12.0, rep.lambda,   1e-9, "total lambda"));
    PetscCall(PicurvAssertIntEqual(0, (PetscInt)rep.cclass, "interior cell class"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityQuickInterior()

static PetscErrorCode TestMomentumStabilityQuickInterior ( void  )

static

Definition at line 699 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->central = 0; simCtx->invicid = 1;       /* QUICK branch */
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 3.0, 0.0, 0.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(4.0, rep.lambda_c, 1e-9, "QUICK interior f_c=4/3 (4/3*3=4)"));
    PetscCall(PicurvAssertIntEqual(0, (PetscInt)rep.cclass, "interior cell class"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityQuickBoundary()

static PetscErrorCode TestMomentumStabilityQuickBoundary ( void  )

static

Definition at line 714 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, INLET));   /* non-periodic -> boundary band active */
    simCtx->central = 0; simCtx->invicid = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 8.0, 0.0, 0.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(20.0, rep.lambda_c, 1e-9, "boundary QUICK f_c=2.5 (2.5*8=20)"));
    PetscCall(PicurvAssertIntEqual(1, (PetscInt)rep.cclass, "physical-boundary cell class"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityViscousLongitudinal()

static PetscErrorCode TestMomentumStabilityViscousLongitudinal ( void  )

static

Definition at line 729 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;          /* nu = 1, no convection */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(24.0, rep.lambda_v, 1e-9, "viscous = 4*6*nu = 24nu (factor-8 longitudinal)"));
    PetscCall(PicurvAssertRealNear(34.0, rep.lambda,   1e-9, "total lambda_t+lambda_v"));
    PetscCall(PicurvAssertIntEqual((PetscInt)MOM_STAB_LIMITER_VISCOSITY, (PetscInt)rep.limiter, "viscosity-limited"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityViscousScalesWithNu()

static PetscErrorCode TestMomentumStabilityViscousScalesWithNu ( void  )

static

Definition at line 744 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 4.0;          /* nu = 0.25 -> 24*0.25 = 6 */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(6.0, rep.lambda_v, 1e-9, "viscous scales with nu=1/Re"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityLESEddyViscosity()

static PetscErrorCode TestMomentumStabilityLESEddyViscosity ( void  )

static

Definition at line 757 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0; simCtx->les = 1;   /* nu_eff = 1 + nu_t */
    PetscCall(MomFillLocalScalar(user->da, user->lNu_t, 1.0)); /* nu_t = 1 -> nu_eff = 2 -> 24*2 = 48 */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(48.0, rep.lambda_v, 1e-9, "LES nu_t adds: 4*6*(1+1) = 48"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityOneSidedViscous()

static PetscErrorCode TestMomentumStabilityOneSidedViscous ( void  )

static

Definition at line 771 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;
    /* Mark one interior cell solid; its fluid cross-neighbours get the one-sided x2 (24 -> 48). */
    PetscCall(MomSetLocalScalarCell(user->da, user->lNvert, 4, 5, 4, 1.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(48.0, rep.lambda_v, 1e-9, "one-sided viscous x2 (24 -> 48)"));
    PetscCall(PicurvAssertIntEqual(1, (PetscInt)rep.one_sided, "one-sided multiplier flagged"));
    PetscCall(PicurvAssertIntEqual(2, (PetscInt)rep.cclass, "IB-adjacent cell class"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityReadOnly()

static PetscErrorCode TestMomentumStabilityReadOnly ( void  )

static

Definition at line 788 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscReal n_ucont_0, n_ucat_0, n_nvert_0, n_ucont_1, n_ucat_1, n_nvert_1;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 2.0;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 1.5, -0.7, 0.3));
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcat,  0.9,  0.4, -0.2));
    PetscCall(VecNorm(user->lUcont, NORM_2, &n_ucont_0));
    PetscCall(VecNorm(user->lUcat,  NORM_2, &n_ucat_0));
    PetscCall(VecNorm(user->lNvert, NORM_2, &n_nvert_0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_D, &rep));
    PetscCall(VecNorm(user->lUcont, NORM_2, &n_ucont_1));
    PetscCall(VecNorm(user->lUcat,  NORM_2, &n_ucat_1));
    PetscCall(VecNorm(user->lNvert, NORM_2, &n_nvert_1));
    PetscCall(PicurvAssertRealNear(n_ucont_0, n_ucont_1, 1e-14, "lUcont unchanged by estimator"));
    PetscCall(PicurvAssertRealNear(n_ucat_0,  n_ucat_1,  1e-14, "lUcat unchanged by estimator"));
    PetscCall(PicurvAssertRealNear(n_nvert_0, n_nvert_1, 1e-14, "lNvert unchanged by estimator"));
    PetscCall(PicurvAssertBool((PetscBool)(rep.lambda > 0.0), "estimate finite and positive"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomSetCartesianMetrics()

static PetscErrorCode MomSetCartesianMetrics ( UserCtx *  user, PetscReal  dx, PetscReal  dy, PetscReal  dz  )

static

Sets uniform Cartesian metrics for spacing (dx,dy,dz) on all cell and face arrays.

Parameters

user	Test user context.
dx,dy,dz	Cartesian cell spacings.

Returns

PetscErrorCode 0 on success.

Definition at line 818 of file test_solver_kernels.c.

{
    const PetscReal sx = dy*dz, sy = dx*dz, sz = dx*dy, aj = 1.0/(dx*dy*dz);
    Vec csiv[] = {user->lCsi, user->lICsi, user->lJCsi, user->lKCsi};
    Vec etav[] = {user->lEta, user->lIEta, user->lJEta, user->lKEta};
    Vec zetv[] = {user->lZet, user->lIZet, user->lJZet, user->lKZet};
    Vec ajv[]  = {user->lAj,  user->lIAj,  user->lJAj,  user->lKAj};
    PetscFunctionBeginUser;
    for (int t = 0; t < 4; ++t) {
        PetscCall(MomFillLocalCmpnts(user->fda, csiv[t], sx, 0, 0));
        PetscCall(MomFillLocalCmpnts(user->fda, etav[t], 0, sy, 0));
        PetscCall(MomFillLocalCmpnts(user->fda, zetv[t], 0, 0, sz));
        PetscCall(MomFillLocalScalar(user->da, ajv[t], aj));
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomFillUcontXRamp()

static PetscErrorCode MomFillUcontXRamp ( DM  fda, Vec  lvec, PetscReal  slope  )

static

Fills lUcont.x as a linear ramp slope*i (nonzero discrete contravariant divergence).

Parameters

fda	Vector DM.
lvec	Local Cmpnts vector.
slope	Ramp slope.

Returns

PetscErrorCode 0 on success.

Definition at line 842 of file test_solver_kernels.c.

{
    Cmpnts ***a;
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(fda, &gxs, &gys, &gzs, &gxm, &gym, &gzm));
    PetscCall(DMDAVecGetArray(fda, lvec, &a));
    for (PetscInt k = gzs; k < gzs+gzm; ++k)
        for (PetscInt j = gys; j < gys+gym; ++j)
            for (PetscInt i = gxs; i < gxs+gxm; ++i) { a[k][j][i].x = slope*i; a[k][j][i].y = 0; a[k][j][i].z = 0; }
    PetscCall(DMDAVecRestoreArray(fda, lvec, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MomFillUcatShearY()

static PetscErrorCode MomFillUcatShearY ( DM  fda, Vec  lvec, PetscReal  gamma  )

static

Fills lUcat as a linear shear u=(gamma*j,0,0) so |grad u|_inf = gamma.

Parameters

fda	Vector DM.
lvec	Local Cmpnts vector.
gamma	Shear rate.

Returns

PetscErrorCode 0 on success.

Definition at line 863 of file test_solver_kernels.c.

{
    Cmpnts ***a;
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(fda, &gxs, &gys, &gzs, &gxm, &gym, &gzm));
    PetscCall(DMDAVecGetArray(fda, lvec, &a));
    for (PetscInt k = gzs; k < gzs+gzm; ++k)
        for (PetscInt j = gys; j < gys+gym; ++j)
            for (PetscInt i = gxs; i < gxs+gxm; ++i) { a[k][j][i].x = gamma*j; a[k][j][i].y = 0; a[k][j][i].z = 0; }
    PetscCall(DMDAVecRestoreArray(fda, lvec, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestMomentumStabilityAnisotropicViscous()

static PetscErrorCode TestMomentumStabilityAnisotropicViscous ( void  )

static

Definition at line 879 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;
    PetscCall(MomSetCartesianMetrics(user, 1.0, 2.0, 4.0));   /* 8*(1 + 1/4 + 1/16) = 10.5 */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(10.5, rep.lambda_v, 1e-9, "anisotropic viscous = 8nu*sum(1/d^2)"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityDirectionalQuick()

static PetscErrorCode TestMomentumStabilityDirectionalQuick ( void  )

static

Definition at line 893 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGridP(&simCtx, &user, 8, PETSC_FALSE, PETSC_TRUE, PETSC_TRUE, INLET)); /* x non-periodic */
    simCtx->central = 0; simCtx->invicid = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 3.0, 3.0, 3.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    /* boundary cell: 0.5*(2.5*6 + (4/3)*6 + (4/3)*6) = 15.5 ; NOT 22.5 (all-2.5). */
    PetscCall(PicurvAssertRealNear(15.5, rep.lambda_c, 1e-9, "only x-direction modified to 2.5"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityOneSidedDirections()

static PetscErrorCode TestMomentumStabilityOneSidedDirections ( void  )

static

Definition at line 908 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    const PetscInt offs[3][3] = {{1,0,0},{0,1,0},{0,0,1}};   /* x,y,z neighbors of (4,4,4) */
    PetscFunctionBeginUser;
    for (int d = 0; d < 3; ++d) {
        PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
        simCtx->invicid = 0; simCtx->ren = 1.0;
        PetscCall(MomSetLocalScalarCell(user->da, user->lNvert,
                                        4+offs[d][0], 4+offs[d][1], 4+offs[d][2], 1.0));
        PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
        PetscCall(PicurvAssertRealNear(48.0, rep.lambda_v, 1e-9, "single-direction one-sided x2"));
        PetscCall(PicurvAssertIntEqual(1, (PetscInt)rep.one_sided, "one-sided flagged"));
        PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    }
    /* multiple solid neighbors still apply x2 only once (48, not 96). */
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;
    PetscCall(MomSetLocalScalarCell(user->da, user->lNvert, 5, 4, 4, 1.0));
    PetscCall(MomSetLocalScalarCell(user->da, user->lNvert, 4, 5, 4, 1.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(48.0, rep.lambda_v, 1e-9, "multiple triggers -> x2 once (not 96)"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityActiveRowMask()

static PetscErrorCode TestMomentumStabilityActiveRowMask ( void  )

static

Definition at line 936 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    /* Positive non-periodic x face (i=mx-2) disables the xi row only; cell stays active. */
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, INLET));
    simCtx->central = 1; simCtx->invicid = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 2.0, 2.0, 2.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertBool((PetscBool)(rep.active_cells > 0), "boundary cells remain active"));
    PetscCall(PicurvAssertBool((PetscBool)(rep.lambda > rep.lambda_t), "tangential rows retained"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    /* All cells solid -> zero active -> estimate falls back to lambda_t only. */
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 1;
    PetscCall(MomFillLocalScalar(user->da, user->lNvert, 1.0));   /* every cell solid */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertIntEqual(0, rep.active_cells, "all solid -> zero active cells"));
    PetscCall(PicurvAssertRealNear(rep.lambda_t, rep.lambda, 1e-12, "empty active set -> lambda_t only"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityValidation()

static PetscErrorCode TestMomentumStabilityValidation ( void  )

static

Definition at line 960 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscErrorCode e1 = 0, e2 = 0;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    e1 = ComputeMomentumStabilityEstimate(user, 0, simCtx->dt, MOM_STAB_CAND_C, &rep);  /* block_number=0 */
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(e1 != 0), "block_number<=0 rejected"));
 
    PetscCall(MomSetLocalScalarCell(user->da, user->lAj, 4, 4, 4, PETSC_INFINITY));     /* bad metric */
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    e2 = ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(e2 != 0), "non-finite inverse-Jacobian rejected"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityCandidateBCDivergence()

static PetscErrorCode TestMomentumStabilityCandidateBCDivergence ( void  )

static

Definition at line 984 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->central = 1; simCtx->invicid = 1;                      /* f_c=1 isolates the div term */
    PetscCall(MomFillUcontXRamp(user->fda, user->lUcont, 2.0));    /* div = slope = 2 */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    /* C - B = 0.5*Aj*|div| = 0.5*2 = 1.0 at the (shared) controlling cell. */
    PetscCall(PicurvAssertRealNear(1.0, rep.lambda_C - rep.lambda_B, 1e-9, "C-B = divergence term"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityCandidateDShear()

static PetscErrorCode TestMomentumStabilityCandidateDShear ( void  )

static

Definition at line 999 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->central = 1; simCtx->invicid = 1;
    PetscCall(MomFillUcatShearY(user->fda, user->lUcat, 2.0));     /* |grad u|_inf = gamma = 2 */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_D, &rep));
    PetscCall(PicurvAssertRealNear(2.0, rep.lambda_D - rep.lambda_C, 1e-9, "D-C = |grad u| on shear"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
 
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->central = 1; simCtx->invicid = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcat, 0.7, -0.3, 0.5));  /* uniform velocity */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_D, &rep));
    PetscCall(PicurvAssertRealNear(rep.lambda_C, rep.lambda_D, 1e-12, "D == C for uniform velocity"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityWallSuppression()

static PetscErrorCode TestMomentumStabilityWallSuppression ( void  )

static

Definition at line 1020 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscFunctionBeginUser;
    /* n=3 -> DMDA mx=4 -> interior cells {1,2}; every interior cell touches 3 wall faces
       (no wall-free interior cell exists, so the global max reflects suppression). */
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 3, WALL));
    simCtx->invicid = 0; simCtx->ren = 1.0; simCtx->les = 1;
    PetscCall(MomFillLocalScalar(user->da, user->lNu_t, 1.0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    /* 3 wall faces (nu_t suppressed, nu_eff=1) + 3 interior faces (nu_eff=2): 4*(3*1+3*2)=36. */
    PetscCall(PicurvAssertRealNear(36.0, rep.lambda_v, 1e-9, "wall-face eddy viscosity suppressed"));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityReadOnlyExact()

static PetscErrorCode TestMomentumStabilityReadOnlyExact ( void  )

static

Definition at line 1037 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    Vec ucont0, ucat0, nvert0, nut0, aj0;
    PetscBool eq;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 2.0; simCtx->les = 1;
    PetscCall(MomFillLocalCmpnts(user->fda, user->lUcont, 1.5, -0.7, 0.3));
    PetscCall(MomFillUcatShearY(user->fda, user->lUcat, 1.1));
    PetscCall(MomFillLocalScalar(user->da, user->lNu_t, 0.4));
    PetscCall(VecDuplicate(user->lUcont, &ucont0)); PetscCall(VecCopy(user->lUcont, ucont0));
    PetscCall(VecDuplicate(user->lUcat,  &ucat0));  PetscCall(VecCopy(user->lUcat,  ucat0));
    PetscCall(VecDuplicate(user->lNvert, &nvert0)); PetscCall(VecCopy(user->lNvert, nvert0));
    PetscCall(VecDuplicate(user->lNu_t,  &nut0));   PetscCall(VecCopy(user->lNu_t,  nut0));
    PetscCall(VecDuplicate(user->lAj,    &aj0));    PetscCall(VecCopy(user->lAj,    aj0));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_D, &rep));
    PetscCall(VecEqual(user->lUcont, ucont0, &eq)); PetscCall(PicurvAssertBool(eq, "lUcont bit-identical"));
    PetscCall(VecEqual(user->lUcat,  ucat0,  &eq)); PetscCall(PicurvAssertBool(eq, "lUcat bit-identical"));
    PetscCall(VecEqual(user->lNvert, nvert0, &eq)); PetscCall(PicurvAssertBool(eq, "lNvert bit-identical"));
    PetscCall(VecEqual(user->lNu_t,  nut0,   &eq)); PetscCall(PicurvAssertBool(eq, "lNu_t bit-identical"));
    PetscCall(VecEqual(user->lAj,    aj0,    &eq)); PetscCall(PicurvAssertBool(eq, "lAj bit-identical"));
    PetscCall(VecDestroy(&ucont0)); PetscCall(VecDestroy(&ucat0)); PetscCall(VecDestroy(&nvert0));
    PetscCall(VecDestroy(&nut0));   PetscCall(VecDestroy(&aj0));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumActiveRowsHelper()

static PetscErrorCode TestMomentumActiveRowsHelper ( void  )

static

Definition at line 1067 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL;
    PetscReal ***nvert;
    const PetscInt mx = 9, my = 9, mz = 9;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));   /* mx=9; all fluid */
    PetscCall(DMDAVecGetArray(user->da, user->lNvert, &nvert));
    /* fully periodic interior cell (4,4,4): all rows active regardless of TwoD value below. */
    PetscCall(PicurvAssertIntEqual(0x7, MomCellActiveRows(nvert,4,4,4,mx,my,mz,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE,0), "all rows active (periodic interior)"));
    PetscCall(PicurvAssertIntEqual(0x6, MomCellActiveRows(nvert,4,4,4,mx,my,mz,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE,1), "TwoD=1 clears xi row"));
    PetscCall(PicurvAssertIntEqual(0x5, MomCellActiveRows(nvert,4,4,4,mx,my,mz,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE,2), "TwoD=2 clears eta row"));
    PetscCall(PicurvAssertIntEqual(0x3, MomCellActiveRows(nvert,4,4,4,mx,my,mz,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE,3), "TwoD=3 clears zeta row"));
    /* positive non-periodic faces clear the corresponding normal rows. */
    PetscCall(PicurvAssertIntEqual(0x6, MomCellActiveRows(nvert,4,4,mx-2,mx,my,mz,PETSC_TRUE,PETSC_FALSE,PETSC_FALSE,0), "pos non-periodic xi face clears xi row"));
    /* a location is skipped only when all three rows are inactive (here: solid cell). */
    PetscCall(MomSetLocalScalarCell(user->da, user->lNvert, 4, 4, 4, 1.0));
    PetscCall(DMDAVecRestoreArray(user->da, user->lNvert, &nvert));
    PetscCall(DMDAVecGetArray(user->da, user->lNvert, &nvert));
    PetscCall(PicurvAssertIntEqual(0, MomCellActiveRows(nvert,4,4,4,mx,my,mz,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE,0), "solid cell -> all rows inactive"));
    PetscCall(DMDAVecRestoreArray(user->da, user->lNvert, &nvert));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumStabilityErrorReentry()

static PetscErrorCode TestMomentumStabilityErrorReentry ( void  )

static

Definition at line 1093 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; MomStabilityReport rep;
    PetscErrorCode e_bad = 0; PetscBool eq;
    Vec aj_good;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 8, PERIODIC));
    simCtx->invicid = 0; simCtx->ren = 1.0;
    PetscCall(VecDuplicate(user->lAj, &aj_good)); PetscCall(VecCopy(user->lAj, aj_good));
 
    PetscCall(MomSetLocalScalarCell(user->da, user->lAj, 4, 4, 4, PETSC_INFINITY));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    e_bad = ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(e_bad != 0), "bad metric triggers error"));
 
    /* Arrays must have been restored: lAj is still writable/usable -> repair it. */
    PetscCall(VecCopy(aj_good, user->lAj));
    PetscCall(VecEqual(user->lAj, aj_good, &eq));
    PetscCall(PicurvAssertBool(eq, "lAj usable after error (arrays were restored)"));
 
    /* Re-entry must now succeed on the repaired state. */
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    PetscCall(PicurvAssertRealNear(24.0, rep.lambda_v, 1e-9, "estimator reusable after repair"));
    PetscCall(VecDestroy(&aj_good));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestMomentumBDFResidualCombination()

static PetscErrorCode TestMomentumBDFResidualCombination ( void  )

static

Definition at line 1124 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL;
    Vec Rhs, ref; PetscReal dt, nrm;
    PetscFunctionBeginUser;
    PetscCall(MomMakeUnitGrid(&simCtx, &user, 6, PERIODIC));
    dt = simCtx->dt;
    PetscCall(VecSet(user->Ucont,     2.0));   /* U^{n}    */
    PetscCall(VecSet(user->Ucont_o,   1.0));   /* U^{n-1}  */
    PetscCall(VecSet(user->Ucont_rm1, 0.5));   /* U^{n-2}  */
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(VecDuplicate(user->Ucont, &ref));
 
    /* ---- BDF2 (step != 1, step != StartStep) ---- */
    simCtx->StartStep = 0; simCtx->step = 5;
    PetscCall(PicurvAssertBool(MomentumUsesBDF2(simCtx), "step5 uses BDF2"));
    {
        const PetscReal a0 = MomentumBDFCoefficient(simCtx);   /* == 1.5 */
        PetscCall(VecSet(Rhs, 3.0));                            /* mock R_spatial */
        PetscCall(VecAXPY(Rhs, -a0/dt,  user->Ucont));
        PetscCall(VecAXPY(Rhs, +2.0/dt, user->Ucont_o));
        PetscCall(VecAXPY(Rhs, -0.5/dt, user->Ucont_rm1));
        /* independent reference: 3 + (-1.5*2 + 2*1 - 0.5*0.5)/dt = 3 + (-3+2-0.25)/0.1 = 3 - 12.5 = -9.5 */
        PetscCall(VecSet(ref, -9.5));
        PetscCall(VecAXPY(ref, -1.0, Rhs)); PetscCall(VecNorm(ref, NORM_INFINITY, &nrm));
        PetscCall(PicurvAssertRealNear(0.0, nrm, 1e-9, "BDF2 temporal vector combination"));
    }
    /* ---- BDF1 (step == 1) ---- */
    simCtx->StartStep = 0; simCtx->step = 1;
    PetscCall(PicurvAssertBool((PetscBool)(!MomentumUsesBDF2(simCtx)), "step1 uses BDF1"));
    {
        const PetscReal a0 = MomentumBDFCoefficient(simCtx);   /* == 1.0 */
        PetscCall(VecSet(Rhs, 3.0));
        PetscCall(VecAXPY(Rhs, -a0/dt,  user->Ucont));
        PetscCall(VecAXPY(Rhs, +1.0/dt, user->Ucont_o));
        /* reference: 3 + (-1*2 + 1*1)/0.1 = 3 + (-1)/0.1 = 3 - 10 = -7.0 */
        PetscCall(VecSet(ref, -7.0));
        PetscCall(VecAXPY(ref, -1.0, Rhs)); PetscCall(VecNorm(ref, NORM_INFINITY, &nrm));
        PetscCall(PicurvAssertRealNear(0.0, nrm, 1e-9, "BDF1 temporal vector combination"));
    }
    PetscCall(VecDestroy(&Rhs)); PetscCall(VecDestroy(&ref));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

Runs the unit-solver PETSc test binary.

Definition at line 1173 of file test_solver_kernels.c.

{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"les-filter-paths", TestLESTestFilterPaths},
        {"analytical-geometry-selection", TestAnalyticalGeometrySelection},
        {"analytical-scalar-verification-helpers", TestAnalyticalScalarVerificationHelpers},
        {"analytical-solution-engine-dispatch", TestAnalyticalSolutionEngineDispatch},
        {"analytical-solution-engine-taylor-green-samples", TestAnalyticalSolutionEngineTaylorGreenSamples},
        {"analytical-solution-for-particles-dispatch", TestAnalyticalSolutionForParticlesDispatch},
        {"compute-eddy-viscosity-les-deterministic-field", TestComputeEddyViscosityLESDeterministicField},
        {"flow-solver-rejects-unsupported-momentum-solver-type", TestFlowSolverRejectsUnsupportedMomentumSolverType},
        {"driven-channel-flow-source", TestDrivenChannelFlowSource},
        {"momentum-bdf-coefficient", TestMomentumBDFCoefficient},
        {"momentum-stability-centered-convection", TestMomentumStabilityCenteredConvection},
        {"momentum-stability-quick-interior", TestMomentumStabilityQuickInterior},
        {"momentum-stability-quick-boundary", TestMomentumStabilityQuickBoundary},
        {"momentum-stability-viscous-cartesian-factor8-estimate", TestMomentumStabilityViscousLongitudinal},
        {"momentum-stability-viscous-scales-with-nu", TestMomentumStabilityViscousScalesWithNu},
        {"momentum-stability-les-eddy-viscosity", TestMomentumStabilityLESEddyViscosity},
        {"momentum-stability-one-sided-viscous", TestMomentumStabilityOneSidedViscous},
        {"momentum-stability-read-only", TestMomentumStabilityReadOnly},
        {"momentum-stability-anisotropic-viscous", TestMomentumStabilityAnisotropicViscous},
        {"momentum-stability-directional-quick", TestMomentumStabilityDirectionalQuick},
        {"momentum-stability-one-sided-directions", TestMomentumStabilityOneSidedDirections},
        {"momentum-stability-active-row-mask", TestMomentumStabilityActiveRowMask},
        {"momentum-stability-validation", TestMomentumStabilityValidation},
        {"momentum-stability-candidate-bc-divergence", TestMomentumStabilityCandidateBCDivergence},
        {"momentum-stability-candidate-d-shear", TestMomentumStabilityCandidateDShear},
        {"momentum-stability-wall-suppression", TestMomentumStabilityWallSuppression},
        {"momentum-stability-read-only-exact", TestMomentumStabilityReadOnlyExact},
        {"momentum-active-rows-helper", TestMomentumActiveRowsHelper},
        {"momentum-stability-error-reentry", TestMomentumStabilityErrorReentry},
        {"momentum-bdf-residual-combination", TestMomentumBDFResidualCombination},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv solver utility tests");
    if (ierr) {
        return (int)ierr;
    }
 
    ierr = PicurvRunTests("unit-solver", cases, sizeof(cases) / sizeof(cases[0]));
    if (ierr) {
        PetscFinalize();
        return (int)ierr;
    }
 
    ierr = PetscFinalize();
    return (int)ierr;
}

Here is the call graph for this function: