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.

{
    double values[3][3][3];
    double weights[3][3][3];
 
    PetscFunctionBeginUser;
    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;
            }
        }
    }
 
    PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(LES_TEST_FILTER_SIMPSON_IK, values, weights), 1.0e-12,
                                   "Simpson-rule filter should preserve a constant field"));
 
    PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(LES_TEST_FILTER_VOLUME_WEIGHTED_BOX, 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(LES_TEST_FILTER_VOLUME_WEIGHTED_BOX, 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 54 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 105 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, PARTICLE_FIELD_ID_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 179 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 258 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 329 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 407 of file test_solver_kernels.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Cmpnts ***ucat = NULL;
    PetscReal ***nu_t = NULL;
    /* CS stores the coefficient C that multiplies Delta^2 |S|, so C = 0.25 is the
       Cs = 0.5 of the classical form. With Aj = 1 the filter width is 1, and the
       linear field below gives |S| = sqrt(2). */
    const PetscReal coefficient   = 0.25;
    const PetscReal expected_nu_t = 0.25 * PetscSqrtReal(2.0);
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
    simCtx->les = DYNAMIC_SMAGORINSKY;
    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, coefficient));
    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 465 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 488 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 552 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 569 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 590 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 612 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 628 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 at restart setup step"));
    simCtx->step = 6;
    PetscCall(PicurvAssertRealNear(1.0, MomentumBDFCoefficient(simCtx), 1e-12, "a0 BDF1 on first restart solve (ti==StartStep+1)"));
    simCtx->step = 7;
    PetscCall(PicurvAssertRealNear(1.5, MomentumBDFCoefficient(simCtx), 1e-12, "a0 BDF2 after first restart solve"));
    simCtx->restartHistoryAvailable = PETSC_TRUE;
    simCtx->step = 5;
    PetscCall(PicurvAssertRealNear(1.5, MomentumBDFCoefficient(simCtx), 1e-12,
                                   "restored history enables BDF2 at restart setup"));
    simCtx->step = 6;
    PetscCall(PicurvAssertRealNear(1.5, MomentumBDFCoefficient(simCtx), 1e-12,
                                   "restored history enables BDF2 on first restart solve"));
    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 659 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)->central = 0; (*simCtx)->invicid = 0; (*simCtx)->block_number = 1;
    (*simCtx)->TwoD = 0; (*simCtx)->les_gradient_model = 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. */
    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 685 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 694 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 712 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 727 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 742 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 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 = 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 770 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 784 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 801 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 831 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 855 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 876 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 892 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 906 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 921 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 949 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 973 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 997 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 1012 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 1033 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 1050 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 1080 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 1106 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 1137 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 1186 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: