test_poisson_rhs.c File Reference

C unit tests for Poisson, RHS, body-force, and diffusivity helpers. More...

#include "test_support.h"
#include "poisson.h"
#include "rhs.h"
#include "setup.h"
#include "verification_sources.h"

Include dependency graph for test_poisson_rhs.c:

Go to the source code of this file.

Functions
static PetscErrorCode	EnsurePoissonAndRhsVectors (UserCtx *user)
	Allocates Poisson/RHS support vectors required by the tests.
static PetscErrorCode	TestUpdatePressureAddsPhi (void)
	Tests that pressure updates add the correction potential.
static PetscErrorCode	TestPoissonRHSZeroDivergence (void)
	Tests that the Poisson RHS is zero for zero-divergence velocity.
static PetscErrorCode	TestComputeBodyForcesDispatcher (void)
	Tests the body-force dispatcher across supported source terms.
static PetscErrorCode	TestComputeEulerianDiffusivityMolecularOnly (void)
	Tests Eulerian diffusivity for the molecular-only configuration.
static PetscErrorCode	TestConvectionZeroField (void)
	Tests that convection vanishes for a quiescent field.
static PetscErrorCode	TestViscousUniformField (void)
	Tests that viscous terms vanish for a uniform field.
static PetscErrorCode	ScaleMinimalMetricsToSpacing (UserCtx *user, PetscReal h)
	Rescales the minimal fixture's unit metrics to cubic cells of edge h.
static PetscErrorCode	ClarkToMolecularRatioPerX (PetscReal h, PetscInt probe, PetscReal *ratio_per_x)
	Returns the Clark contribution over the molecular viscous one, per unit x, at one cell.
static PetscErrorCode	TestClarkGradientModelScalesWithFilterWidthSquared (void)
	Tests that the Clark gradient model scales as the square of the filter width.
static PetscErrorCode	TestComputeRHSZeroFieldNoForcing (void)
	Tests that the full RHS remains zero without forcing on a quiescent field.
static PetscErrorCode	TestComputeEulerianDiffusivityGradientConstantField (void)
	Tests diffusivity-gradient computation on a constant field.
static PetscErrorCode	TestComputeEulerianDiffusivityVerificationLinearX (void)
	Tests verification-driven linear diffusivity override and its gradient.
static PetscErrorCode	TestAssemblePoissonOperatorPopulatesRows (void)
	Tests that Poisson matrix assembly produces a populated operator on a tiny Cartesian grid.
static PetscErrorCode	RandomizeMinimalFaceMetrics (UserCtx *user)
	Fills every face metric and Jacobian of the minimal fixture with reproducible random values, so each face gradient carries all of its cross terms.
static PetscErrorCode	TestOperatorAndProjectionShareOneFaceGradient (void)
	Tests that the operator is exactly the divergence of the projection's gradient.
static PetscErrorCode	TestProjectionZeroPhiLeavesVelocityUnchanged (void)
	Tests that projection leaves a zero pressure-correction field unchanged.
static PetscErrorCode	TestProjectionLinearPhiCorrectsVelocity (void)
	Tests that projection applies the expected x-direction correction for a linear pressure field.
static PetscErrorCode	PerturbInteriorFaceFluxes (UserCtx *user, PetscReal amplitude)
	Adds a smooth divergent perturbation to every interior face flux.
static PetscErrorCode	TestPoissonSolverMultigridProjectsToDivergenceFree (void)
	Tests that the production multigrid Poisson solve and projection remove a divergence.
static PetscErrorCode	CountFileOccurrences (const char *path, const char *needle, PetscInt *count)
	Counts the occurrences of needle in a text file; a missing file counts zero.
static PetscErrorCode	TestPoissonSolverMultigridReusesItsSolver (void)
	Tests that the multigrid solver is built once and reused on later steps.
static PetscErrorCode	TestPoissonSolverMultigridHonoursDistinctSweeps (void)
	Tests that unequal pre- and post-smoothing counts give two distinct smoothers.
static PetscErrorCode	TestPoissonNullSpaceRemovesTheInteriorMean (void)
	Tests the null space attached to the operator.
static PetscErrorCode	TestPoissonSolverMultigridRefusesAnOvercoarsenedHierarchy (void)
	Tests that an over-coarsened hierarchy stops the run at the Poisson solve.
int	main (int argc, char **argv)
	Runs the unit-poisson-rhs PETSc test binary.

Detailed Description

C unit tests for Poisson, RHS, body-force, and diffusivity helpers.

Definition in file test_poisson_rhs.c.

Function Documentation

EnsurePoissonAndRhsVectors()

static PetscErrorCode EnsurePoissonAndRhsVectors ( UserCtx *  user )

static

Allocates Poisson/RHS support vectors required by the tests.

Definition at line 16 of file test_poisson_rhs.c.

{
    PetscFunctionBeginUser;
    if (!user->Phi) PetscCall(DMCreateGlobalVector(user->da, &user->Phi));
    if (!user->lPhi) PetscCall(DMCreateLocalVector(user->da, &user->lPhi));
    if (!user->Aj) PetscCall(DMCreateGlobalVector(user->da, &user->Aj));
    if (!user->lAj) PetscCall(DMCreateLocalVector(user->da, &user->lAj));
    if (!user->Diffusivity) PetscCall(DMCreateGlobalVector(user->da, &user->Diffusivity));
    if (!user->lDiffusivity) PetscCall(DMCreateLocalVector(user->da, &user->lDiffusivity));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestUpdatePressureAddsPhi()

static PetscErrorCode TestUpdatePressureAddsPhi ( void  )

static

Tests that pressure updates add the correction potential.

Definition at line 31 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(EnsurePoissonAndRhsVectors(user));
    PetscCall(VecSet(user->P, 1.0));
    PetscCall(VecSet(user->Phi, 0.25));
 
    PetscCall(UpdatePressure(user));
    PetscCall(PicurvAssertVecConstant(user->P, 1.25, 1.0e-12, "UpdatePressure should add Phi into P"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestPoissonRHSZeroDivergence()

static PetscErrorCode TestPoissonRHSZeroDivergence ( void  )

static

Tests that the Poisson RHS is zero for zero-divergence velocity.

Definition at line 52 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec B = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(EnsurePoissonAndRhsVectors(user));
 
    simCtx->dt = 0.5;
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Aj, INSERT_VALUES, user->lAj));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Aj, INSERT_VALUES, user->lAj));
 
    PetscCall(VecDuplicate(user->P, &B));
    PetscCall(ComputePoissonRHS(user, B));
    PetscCall(PicurvAssertVecConstant(B, 0.0, 1.0e-12, "zero velocity divergence should produce zero Poisson RHS"));
    PetscCall(PicurvAssertRealNear(0.0, simCtx->poissonSourceImbalance, 1.0e-12, "global Poisson RHS sum should be zero"));
 
    PetscCall(VecDestroy(&B));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestComputeBodyForcesDispatcher()

static PetscErrorCode TestComputeBodyForcesDispatcher ( void  )

static

Tests the body-force dispatcher across supported source terms.

Definition at line 86 of file test_poisson_rhs.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(ComputeBodyForces(user, rct));
    PetscCall(DMDAVecGetArrayRead(user->fda, rct, &rct_arr));
    PetscCall(PicurvAssertRealNear(1.5, simCtx->drivingForceMagnitude, 1.0e-12,
                                   "ComputeBodyForces should update the driven-flow controller magnitude"));
    PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].y, 1.0e-12, "ComputeBodyForces should keep y unchanged"));
    PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].z, 1.0e-12, "ComputeBodyForces should keep z 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:

TestComputeEulerianDiffusivityMolecularOnly()

static PetscErrorCode TestComputeEulerianDiffusivityMolecularOnly ( void  )

static

Tests Eulerian diffusivity for the molecular-only configuration.

Definition at line 140 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscReal ***diff_arr = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(EnsurePoissonAndRhsVectors(user));
 
    simCtx->ren = 2.0;
    simCtx->schmidt_number = 4.0;
    simCtx->les = PETSC_FALSE;
    PetscCall(VecZeroEntries(user->Diffusivity));
 
    PetscCall(ComputeEulerianDiffusivity(user));
    PetscCall(DMDAVecGetArrayRead(user->da, user->Diffusivity, &diff_arr));
    PetscCall(PicurvAssertRealNear(0.125, diff_arr[1][1][1], 1.0e-12, "molecular diffusivity interior value"));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Diffusivity, &diff_arr));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestConvectionZeroField()

static PetscErrorCode TestConvectionZeroField ( void  )

static

Tests that convection vanishes for a quiescent field.

Definition at line 167 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec conv = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecDuplicate(user->lUcont, &conv));
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Ucat, 0.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
 
    PetscCall(Convection(user, user->lUcont, user->lUcat, conv));
    PetscCall(PicurvAssertVecConstant(conv, 0.0, 1.0e-12, "Convection should vanish for a quiescent field"));
 
    PetscCall(VecDestroy(&conv));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestViscousUniformField()

static PetscErrorCode TestViscousUniformField ( void  )

static

Tests that viscous terms vanish for a uniform field.

Definition at line 194 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec visc = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecDuplicate(user->lUcont, &visc));
    PetscCall(VecSet(user->Ucont, 1.0));
    PetscCall(VecSet(user->Ucat, 1.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
 
    PetscCall(Viscous(user, user->lUcont, user->lUcat, visc));
    PetscCall(PicurvAssertVecConstant(visc, 0.0, 1.0e-12, "Viscous should vanish for a uniform field"));
 
    PetscCall(VecDestroy(&visc));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ScaleMinimalMetricsToSpacing()

static PetscErrorCode ScaleMinimalMetricsToSpacing ( UserCtx *  user, PetscReal  h  )

static

Rescales the minimal fixture's unit metrics to cubic cells of edge h.

Definition at line 220 of file test_poisson_rhs.c.

{
    Vec area_vectors[12] = {user->Csi, user->Eta, user->Zet, user->ICsi, user->IEta, user->IZet,
                            user->JCsi, user->JEta, user->JZet, user->KCsi, user->KEta, user->KZet};
    Vec local_area[12]   = {user->lCsi, user->lEta, user->lZet, user->lICsi, user->lIEta, user->lIZet,
                            user->lJCsi, user->lJEta, user->lJZet, user->lKCsi, user->lKEta, user->lKZet};
    Vec jacobians[4]     = {user->Aj, user->IAj, user->JAj, user->KAj};
    Vec local_jac[4]     = {user->lAj, user->lIAj, user->lJAj, user->lKAj};
 
    PetscFunctionBeginUser;
    for (PetscInt v = 0; v < 12; ++v) {
        PetscCall(VecScale(area_vectors[v], h * h));   /* a face of a cube of edge h */
        PetscCall(DMGlobalToLocalBegin(user->fda, area_vectors[v], INSERT_VALUES, local_area[v]));
        PetscCall(DMGlobalToLocalEnd(user->fda, area_vectors[v], INSERT_VALUES, local_area[v]));
    }
    for (PetscInt v = 0; v < 4; ++v) {
        PetscCall(VecSet(jacobians[v], 1.0 / (h * h * h)));
        PetscCall(DMGlobalToLocalBegin(user->da, jacobians[v], INSERT_VALUES, local_jac[v]));
        PetscCall(DMGlobalToLocalEnd(user->da, jacobians[v], INSERT_VALUES, local_jac[v]));
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ClarkToMolecularRatioPerX()

static PetscErrorCode ClarkToMolecularRatioPerX ( PetscReal  h, PetscInt  probe, PetscReal *  ratio_per_x  )

static

Returns the Clark contribution over the molecular viscous one, per unit x, at one cell.

The field is u_x = x^2 on cubic cells of edge h, so du/dx = 2x and the exact gradient-model stress tau_11 = (h^2/12)(du/dx)^2 varies along x. Both the model's and viscosity's contributions leave Viscous() through the same flux differences, so their ratio does not depend on how Viscous() normalizes its output.

Definition at line 252 of file test_poisson_rhs.c.

{
    SimCtx   *simCtx = NULL;
    UserCtx  *user = NULL;
    Vec       molecular = NULL, with_clark = NULL;
    Cmpnts ***ucat = NULL, ***visc0 = NULL, ***visc1 = NULL;
    const PetscInt n = 10, mid = n / 2;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, n, n, n));
    PetscCall(ScaleMinimalMetricsToSpacing(user, h));
    simCtx->ren = 1.0;
    simCtx->les = NO_LES_MODEL;
 
    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) {
        const PetscReal x = (i - 0.5) * h;   /* cell centre */
        ucat[k][j][i].x = x * x; 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(VecSet(user->lUcont, 0.0));
 
    PetscCall(VecDuplicate(user->lUcont, &molecular));
    PetscCall(VecDuplicate(user->lUcont, &with_clark));
    simCtx->les_gradient_model = 0;
    PetscCall(Viscous(user, user->lUcont, user->lUcat, molecular));
    simCtx->les_gradient_model = 1;
    PetscCall(Viscous(user, user->lUcont, user->lUcat, with_clark));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, molecular, &visc0));
    PetscCall(DMDAVecGetArrayRead(user->fda, with_clark, &visc1));
    {
        const PetscReal v_molecular = visc0[mid][mid][probe].x;
        const PetscReal v_clark     = visc1[mid][mid][probe].x - v_molecular;
        PetscCheck(PetscAbsReal(v_molecular) > 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB,
                   "The molecular viscous term vanished at the probe; the fixture is wrong.");
        *ratio_per_x = (v_clark / v_molecular) / ((probe - 0.5) * h);
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, with_clark, &visc1));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, molecular, &visc0));
    PetscCall(VecDestroy(&with_clark));
    PetscCall(VecDestroy(&molecular));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestClarkGradientModelScalesWithFilterWidthSquared()

static PetscErrorCode TestClarkGradientModelScalesWithFilterWidthSquared ( void  )

static

Tests that the Clark gradient model scales as the square of the filter width.

The gradient-model stress is tau_ij = (Delta_k^2/12) du_i/dx_k du_j/dx_k, so for a fixed physical field its force grows as Delta^2 relative to the molecular viscous force. Halving the spacing must therefore divide the Clark-to-molecular ratio at a given x by four. Fixtures with unit cells hide any error in the power of Delta, which is why the spacing is varied here.

Definition at line 311 of file test_poisson_rhs.c.

{
    PetscReal coarse = 0.0, fine = 0.0;
 
    PetscFunctionBeginUser;
    /* Probe cells chosen so the two spacings give comparable x without being on a
       boundary: cell 5 at h = 0.2 and cell 5 at h = 0.1. The ratio is per unit x. */
    PetscCall(ClarkToMolecularRatioPerX(0.2, 5, &coarse));
    PetscCall(ClarkToMolecularRatioPerX(0.1, 5, &fine));
    PetscCall(PicurvAssertBool((PetscBool)(coarse < 0.0 && fine < 0.0),
                               "the gradient model removes momentum where du/dx grows along x"));
    PetscCall(PicurvAssertRealNear(4.0, coarse / fine, 1.0e-8,
                                   "halving the filter width divides the gradient-model force by four"));
    /* Exact value: the molecular flux carries both halves of the symmetric stress,
       2 nu du/dx per unit area, so its difference across a cell is 4 nu h^3 for
       u = x^2; the model's is -(2/3) h^5 x. Their ratio per unit x is -h^2/(6 nu),
       which pins the 1/12 coefficient as well as the power of Delta. */
    PetscCall(PicurvAssertRealNear(-0.2 * 0.2 / 6.0, coarse, 1.0e-10,
                                   "gradient-model force matches (Delta^2/12) du/dx du/dx at h = 0.2"));
    PetscCall(PicurvAssertRealNear(-0.1 * 0.1 / 6.0, fine, 1.0e-10,
                                   "gradient-model force matches (Delta^2/12) du/dx du/dx at h = 0.1"));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestComputeRHSZeroFieldNoForcing()

static PetscErrorCode TestComputeRHSZeroFieldNoForcing ( void  )

static

Tests that the full RHS remains zero without forcing on a quiescent field.

Definition at line 339 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Ucat, 0.0));
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    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->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
 
    PetscCall(ComputeRHS(user, user->Rhs));
    PetscCall(PicurvAssertVecConstant(user->Rhs, 0.0, 1.0e-12, "ComputeRHS should remain zero without forcing on a quiescent field"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestComputeEulerianDiffusivityGradientConstantField()

static PetscErrorCode TestComputeEulerianDiffusivityGradientConstantField ( void  )

static

Tests diffusivity-gradient computation on a constant field.

Definition at line 366 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->Diffusivity, 0.75));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Diffusivity, INSERT_VALUES, user->lDiffusivity));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Diffusivity, INSERT_VALUES, user->lDiffusivity));
 
    PetscCall(ComputeEulerianDiffusivityGradient(user));
    PetscCall(PicurvAssertVecConstant(user->DiffusivityGradient, 0.0, 1.0e-12, "constant diffusivity should yield zero diffusivity gradient"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestComputeEulerianDiffusivityVerificationLinearX()

static PetscErrorCode TestComputeEulerianDiffusivityVerificationLinearX ( void  )

static

Tests verification-driven linear diffusivity override and its gradient.

Definition at line 387 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscReal ***diff_arr = NULL;
    Cmpnts ***grad_arr = NULL;
    Cmpnts ***cent = NULL;
    const PetscReal gamma0 = 0.5;
    const PetscReal slope_x = 0.25;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PetscStrncpy(simCtx->eulerianSource, "analytical", sizeof(simCtx->eulerianSource)));
    simCtx->verificationDiffusivity.enabled = PETSC_TRUE;
    PetscCall(PetscStrncpy(simCtx->verificationDiffusivity.mode, "analytical", sizeof(simCtx->verificationDiffusivity.mode)));
    PetscCall(PetscStrncpy(simCtx->verificationDiffusivity.profile, "LINEAR_X", sizeof(simCtx->verificationDiffusivity.profile)));
    simCtx->verificationDiffusivity.gamma0 = gamma0;
    simCtx->verificationDiffusivity.slope_x = slope_x;
 
    PetscCall(PicurvAssertBool(VerificationDiffusivityOverrideActive(simCtx),
                               "verification diffusivity override should report as active"));
    PetscCall(ComputeEulerianDiffusivity(user));
    PetscCall(ComputeEulerianDiffusivityGradient(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->Diffusivity, &diff_arr));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Cent, &cent));
    PetscCall(PicurvAssertRealNear(gamma0 + slope_x * cent[2][2][2].x, diff_arr[2][2][2], 1.0e-12,
                                   "verification override should populate the linear diffusivity field"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Cent, &cent));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Diffusivity, &diff_arr));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->DiffusivityGradient, &grad_arr));
    PetscCall(PicurvAssertRealNear(slope_x / user->IM, grad_arr[2][2][2].x, 1.0e-12,
                                   "linear verification diffusivity should yield the current finite-difference x gradient"));
    PetscCall(PicurvAssertRealNear(0.0, grad_arr[2][2][2].y, 1.0e-12,
                                   "linear verification diffusivity should keep y gradient zero"));
    PetscCall(PicurvAssertRealNear(0.0, grad_arr[2][2][2].z, 1.0e-12,
                                   "linear verification diffusivity should keep z gradient zero"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->DiffusivityGradient, &grad_arr));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestAssemblePoissonOperatorPopulatesRows()

static PetscErrorCode TestAssemblePoissonOperatorPopulatesRows ( void  )

static

Tests that Poisson matrix assembly produces a populated operator on a tiny Cartesian grid.

Definition at line 433 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscInt rows = 0, cols = 0;
    PetscInt ncols = 0;
    const PetscInt *col_idx = NULL;
    const PetscScalar *values = NULL;
    PetscInt interior_row = 0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    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(AssemblePoissonOperator(user));
    PetscCall(PicurvAssertBool((PetscBool)(user->A != NULL), "AssemblePoissonOperator should allocate the Poisson operator"));
 
    PetscCall(MatGetSize(user->A, &rows, &cols));
    PetscCall(PicurvAssertIntEqual(user->info.mx * user->info.my * user->info.mz, rows, "Poisson operator row count should match the DA node count"));
    PetscCall(PicurvAssertIntEqual(rows, cols, "Poisson operator should be square"));
 
    interior_row = (2 * user->info.my + 2) * user->info.mx + 2;
    PetscCall(MatGetRow(user->A, interior_row, &ncols, &col_idx, &values));
    PetscCall(PicurvAssertBool((PetscBool)(ncols > 1), "Interior Poisson row should couple to neighboring nodes"));
    PetscCall(PicurvAssertBool((PetscBool)(values != NULL), "Interior Poisson row should expose non-null coefficients"));
    PetscCall(MatRestoreRow(user->A, interior_row, &ncols, &col_idx, &values));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RandomizeMinimalFaceMetrics()

static PetscErrorCode RandomizeMinimalFaceMetrics ( UserCtx *  user )

static

Fills every face metric and Jacobian of the minimal fixture with reproducible random values, so each face gradient carries all of its cross terms.

Definition at line 469 of file test_poisson_rhs.c.

{
    Vec         face_metrics[9] = {user->ICsi, user->IEta, user->IZet, user->JCsi, user->JEta,
                                   user->JZet, user->KCsi, user->KEta, user->KZet};
    Vec         local_metrics[9] = {user->lICsi, user->lIEta, user->lIZet, user->lJCsi, user->lJEta,
                                    user->lJZet, user->lKCsi, user->lKEta, user->lKZet};
    Vec         jacobians[4] = {user->Aj, user->IAj, user->JAj, user->KAj};
    Vec         local_jac[4] = {user->lAj, user->lIAj, user->lJAj, user->lKAj};
    PetscRandom rnd;
 
    PetscFunctionBeginUser;
    PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rnd));
    PetscCall(PetscRandomSetInterval(rnd, 0.5, 1.5));
    PetscCall(PetscRandomSetSeed(rnd, 2026));
    PetscCall(PetscRandomSeed(rnd));
    for (PetscInt v = 0; v < 9; ++v) {
        PetscCall(VecSetRandom(face_metrics[v], rnd));
        PetscCall(DMGlobalToLocalBegin(user->fda, face_metrics[v], INSERT_VALUES, local_metrics[v]));
        PetscCall(DMGlobalToLocalEnd(user->fda, face_metrics[v], INSERT_VALUES, local_metrics[v]));
    }
    for (PetscInt v = 0; v < 4; ++v) {
        PetscCall(VecSetRandom(jacobians[v], rnd));
        PetscCall(DMGlobalToLocalBegin(user->da, jacobians[v], INSERT_VALUES, local_jac[v]));
        PetscCall(DMGlobalToLocalEnd(user->da, jacobians[v], INSERT_VALUES, local_jac[v]));
    }
    PetscCall(PetscRandomDestroy(&rnd));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestOperatorAndProjectionShareOneFaceGradient()

static PetscErrorCode TestOperatorAndProjectionShareOneFaceGradient ( void  )

static

Tests that the operator is exactly the divergence of the projection's gradient.

With random non-orthogonal metrics and one solid cell, projecting a zero flux with a random Phi and forming the right-hand side of the result must reproduce -A Phi on every fluid row. The identity holds for any metric values only if the operator and the projection use the same face gradient, including its one-sided forms beside walls and solid cells.

Definition at line 506 of file test_poisson_rhs.c.

{
    SimCtx            *simCtx = NULL;
    UserCtx           *user = NULL;
    Vec                B = NULL, APhi = NULL;
    PetscRandom        rnd;
    PetscReal       ***nvert;
    const PetscReal ***b, ***aphi;
    PetscReal          worst = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 7, 7, 7));
    PetscCall(EnsurePoissonAndRhsVectors(user));
    simCtx->dt = 0.3;
    PetscCall(RandomizeMinimalFaceMetrics(user));
 
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(DMDAVecGetArray(user->da, user->Nvert, &nvert));
    if (user->info.xs <= 3 && 3 < user->info.xs + user->info.xm &&
        user->info.ys <= 3 && 3 < user->info.ys + user->info.ym &&
        user->info.zs <= 3 && 3 < user->info.zs + user->info.zm) nvert[3][3][3] = 1.0;
    PetscCall(DMDAVecRestoreArray(user->da, user->Nvert, &nvert));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
 
    PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rnd));
    PetscCall(PetscRandomSetSeed(rnd, 7));
    PetscCall(PetscRandomSeed(rnd));
    PetscCall(VecSetRandom(user->Phi, rnd));
    PetscCall(PetscRandomDestroy(&rnd));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
 
    PetscCall(AssemblePoissonOperator(user));
    PetscCall(ProjectVelocity(user));
    PetscCall(VecDuplicate(user->P, &B));
    PetscCall(ComputePoissonRHS(user, B));
    PetscCall(VecDuplicate(user->Phi, &APhi));
    PetscCall(MatMult(user->A, user->Phi, APhi));
 
    PetscCall(DMDAVecGetArrayRead(user->da, B, &b));
    PetscCall(DMDAVecGetArrayRead(user->da, APhi, &aphi));
    PetscCall(DMDAVecGetArray(user->da, user->Nvert, &nvert));
    for (PetscInt k = PetscMax(user->info.zs, 1); k < PetscMin(user->info.zs + user->info.zm, user->info.mz - 1); k++) {
        for (PetscInt j = PetscMax(user->info.ys, 1); j < PetscMin(user->info.ys + user->info.ym, user->info.my - 1); j++) {
            for (PetscInt i = PetscMax(user->info.xs, 1); i < PetscMin(user->info.xs + user->info.xm, user->info.mx - 1); i++) {
                if (nvert[k][j][i] > 0.1) continue;
                worst = PetscMax(worst, PetscAbsReal(b[k][j][i] + aphi[k][j][i]) / (1.0 + PetscAbsReal(aphi[k][j][i])));
            }
        }
    }
    PetscCall(DMDAVecRestoreArray(user->da, user->Nvert, &nvert));
    PetscCall(DMDAVecRestoreArrayRead(user->da, APhi, &aphi));
    PetscCall(DMDAVecRestoreArrayRead(user->da, B, &b));
    PetscCall(PicurvAssertRealNear(0.0, worst, 1.0e-12,
                                   "the right-hand side of the projected flux should equal -A Phi on every fluid row"));
 
    PetscCall(VecDestroy(&APhi));
    PetscCall(VecDestroy(&B));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestProjectionZeroPhiLeavesVelocityUnchanged()

static PetscErrorCode TestProjectionZeroPhiLeavesVelocityUnchanged ( void  )

static

Tests that projection leaves a zero pressure-correction field unchanged.

Definition at line 575 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Phi, 0.0));
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
 
    PetscCall(ProjectVelocity(user));
    PetscCall(PicurvAssertVecConstant(user->Ucont, 0.0, 1.0e-12, "ProjectVelocity should leave a zero-velocity field unchanged when Phi is zero"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestProjectionLinearPhiCorrectsVelocity()

static PetscErrorCode TestProjectionLinearPhiCorrectsVelocity ( void  )

static

Tests that projection applies the expected x-direction correction for a linear pressure field.

Definition at line 601 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscReal ***phi = NULL;
    Cmpnts ***ucont = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->dt = COEF_TIME_ACCURACY;
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(DMDAVecGetArray(user->da, user->Phi, &phi));
    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) {
                phi[k][j][i] = (PetscReal)i;
            }
        }
    }
    PetscCall(DMDAVecRestoreArray(user->da, user->Phi, &phi));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Phi, INSERT_VALUES, user->lPhi));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
 
    PetscCall(ProjectVelocity(user));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(PicurvAssertRealNear(-1.0, ucont[2][2][2].x, 1.0e-10, "ProjectVelocity should subtract the x pressure gradient under identity metrics"));
    PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].y, 1.0e-10, "ProjectVelocity should leave the y component unchanged for an x-only gradient"));
    PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].z, 1.0e-10, "ProjectVelocity should leave the z component unchanged for an x-only gradient"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PerturbInteriorFaceFluxes()

static PetscErrorCode PerturbInteriorFaceFluxes ( UserCtx *  user, PetscReal  amplitude  )

static

Adds a smooth divergent perturbation to every interior face flux.

Boundary faces are left alone, so the inflow still balances the outflow and the Neumann pressure problem stays solvable.

Definition at line 645 of file test_poisson_rhs.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***ucont;
 
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArray(user->fda, user->Ucont, &ucont));
    for (PetscInt k = PetscMax(info.zs, 1); k < PetscMin(info.zs + info.zm, info.mz - 1); k++) {
        for (PetscInt j = PetscMax(info.ys, 1); j < PetscMin(info.ys + info.ym, info.my - 1); j++) {
            for (PetscInt i = PetscMax(info.xs, 1); i < PetscMin(info.xs + info.xm, info.mx - 1); i++) {
                if (i <= info.mx - 3) ucont[k][j][i].x += amplitude * PetscSinReal(1.3 * i + 0.7 * j + 0.4 * k);
                if (j <= info.my - 3) ucont[k][j][i].y += amplitude * PetscCosReal(0.9 * i + 1.1 * j + 0.3 * k);
                if (k <= info.mz - 3) ucont[k][j][i].z += amplitude * PetscSinReal(0.5 * i + 0.8 * j + 1.7 * k);
            }
        }
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucont, &ucont));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCONT));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestPoissonSolverMultigridProjectsToDivergenceFree()

static PetscErrorCode TestPoissonSolverMultigridProjectsToDivergenceFree ( void  )

static

Tests that the production multigrid Poisson solve and projection remove a divergence.

Runs PoissonSolver_Multigrid through the real setup path with a three-level hierarchy, then the same UpdatePressure/ProjectVelocity pair the time loop uses, and requires the projected field to be divergence-free to the solve tolerance.

Definition at line 672 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
    PetscReal before;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(
        NULL, PETSC_FALSE,
        "-im 17\n-jm 17\n-km 17\n-mg_level 3\n"
        "-ps_ksp_rtol 1.0e-12\n-ps_ksp_atol 1.0e-14\n-ps_ksp_max_it 200\n",
        &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(PicurvAssertIntEqual(3, simCtx->usermg.mglevels, "fixture should build a three-level hierarchy"));
 
    PetscCall(PerturbInteriorFaceFluxes(user, 0.2));
    PetscCall(ComputeDivergence(user));
    before = simCtx->MaxDiv;
    PetscCall(PicurvAssertBool((PetscBool)(before > 1.0e-2), "perturbation should make the field divergent"));
 
    PetscCall(PoissonSolver_Multigrid(&simCtx->usermg));
    PetscCall(UpdatePressure(user));
    PetscCall(ProjectVelocity(user));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCONT));
    PetscCall(ComputeDivergence(user));
    PetscCall(PicurvAssertBool((PetscBool)(simCtx->MaxDiv < 1.0e-9 * before),
                               "multigrid solve plus projection should leave no divergence above the solve tolerance"));
 
    PetscCall(PicurvDestroyRuntimeContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CountFileOccurrences()

static PetscErrorCode CountFileOccurrences ( const char *  path, const char *  needle, PetscInt *  count  )

static

Counts the occurrences of needle in a text file; a missing file counts zero.

Definition at line 708 of file test_poisson_rhs.c.

{
    FILE *file;
    char  line[1024];
 
    PetscFunctionBeginUser;
    *count = 0;
    file = fopen(path, "r");
    if (!file) PetscFunctionReturn(0);
    while (fgets(line, sizeof(line), file)) {
        if (strstr(line, needle)) (*count)++;
    }
    fclose(file);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TestPoissonSolverMultigridReusesItsSolver()

static PetscErrorCode TestPoissonSolverMultigridReusesItsSolver ( void  )

static

Tests that the multigrid solver is built once and reused on later steps.

Two consecutive solves must use the same Krylov solver and operator, both must converge, and the convergence log must hold one header per step.

Definition at line 729 of file test_poisson_rhs.c.

{
    SimCtx            *simCtx = NULL;
    UserCtx           *user = NULL;
    char               tmpdir[PETSC_MAX_PATH_LEN], log_path[PETSC_MAX_PATH_LEN + 64];
    KSP                first_ksp;
    Mat                first_operator;
    KSPConvergedReason reason;
    PetscInt           headers = 0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(
        NULL, PETSC_FALSE,
        "-im 17\n-jm 17\n-km 17\n-mg_level 3\n"
        "-ps_ksp_rtol 1.0e-12\n-ps_ksp_atol 1.0e-14\n-ps_ksp_max_it 200\n",
        &simCtx, &user, tmpdir, sizeof(tmpdir)));
 
    simCtx->step = simCtx->StartStep + 1;
    PetscCall(PerturbInteriorFaceFluxes(user, 0.2));
    PetscCall(PoissonSolver_Multigrid(&simCtx->usermg));
    first_ksp = user->ksp;
    first_operator = user->A;
    PetscCall(PicurvAssertBool((PetscBool)(first_ksp != NULL), "the first solve should build and keep the solver"));
    PetscCall(KSPGetConvergedReason(user->ksp, &reason));
    PetscCall(PicurvAssertBool((PetscBool)(reason > 0), "the first solve should converge"));
 
    simCtx->step++;
    PetscCall(PerturbInteriorFaceFluxes(user, 0.1));
    PetscCall(PoissonSolver_Multigrid(&simCtx->usermg));
    PetscCall(PicurvAssertBool((PetscBool)(user->ksp == first_ksp), "the second solve should reuse the solver"));
    PetscCall(PicurvAssertBool((PetscBool)(user->A == first_operator), "the second solve should reuse the operator"));
    PetscCall(KSPGetConvergedReason(user->ksp, &reason));
    PetscCall(PicurvAssertBool((PetscBool)(reason > 0), "the reused solver should converge"));
 
    PetscCall(PetscSNPrintf(log_path, sizeof(log_path), "%s/Poisson_Solver_Convergence_History_Block_0.log", simCtx->log_dir));
    PetscCall(CountFileOccurrences(log_path, "--- Convergence for Timestep", &headers));
    PetscCall(PicurvAssertIntEqual(2, headers, "the convergence log should hold one header per solve"));
 
    PetscCall(PicurvDestroyRuntimeContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestPoissonSolverMultigridHonoursDistinctSweeps()

static PetscErrorCode TestPoissonSolverMultigridHonoursDistinctSweeps ( void  )

static

Tests that unequal pre- and post-smoothing counts give two distinct smoothers.

With pre_sweeps 1 and post_sweeps 3 every smoothed level must run one sweep before the coarse correction and three after it, the post-smoother must keep the configured smoother type, and the solve must still remove the divergence. With equal counts the shared smoother runs the configured sweeps.

Definition at line 779 of file test_poisson_rhs.c.

{
    const char *const configurations[2] = {
        "-im 17\n-jm 17\n-km 17\n-mg_level 3\n-mg_pre_it 1\n-mg_post_it 3\n"
        "-ps_mg_levels_1_ksp_type richardson\n-ps_mg_levels_2_ksp_type richardson\n"
        "-ps_ksp_rtol 1.0e-12\n-ps_ksp_atol 1.0e-14\n-ps_ksp_max_it 200\n",
        "-im 17\n-jm 17\n-km 17\n-mg_level 3\n-mg_pre_it 2\n-mg_post_it 2\n",
    };
 
    PetscFunctionBeginUser;
    for (PetscInt c = 0; c < 2; c++) {
        SimCtx   *simCtx = NULL;
        UserCtx  *user = NULL;
        char      tmpdir[PETSC_MAX_PATH_LEN];
        PC        pc;
        PetscReal before;
 
        PetscCall(PicurvBuildTinyRuntimeContextWithOptions(NULL, PETSC_FALSE, configurations[c],
                                                           &simCtx, &user, tmpdir, sizeof(tmpdir)));
        PetscCall(PerturbInteriorFaceFluxes(user, 0.2));
        PetscCall(ComputeDivergence(user));
        before = simCtx->MaxDiv;
        PetscCall(PoissonSolver_Multigrid(&simCtx->usermg));
        PetscCall(KSPGetPC(user->ksp, &pc));
 
        for (PetscInt l = 1; l < simCtx->usermg.mglevels; l++) {
            KSP       pre_smoother, post_smoother;
            PetscInt  pre_its, post_its;
            PetscBool is_richardson;
 
            PetscCall(PCMGGetSmootherDown(pc, l, &pre_smoother));
            PetscCall(KSPGetTolerances(pre_smoother, NULL, NULL, NULL, &pre_its));
            if (c == 0) {
                PetscCall(PCMGGetSmootherUp(pc, l, &post_smoother));
                PetscCall(KSPGetTolerances(post_smoother, NULL, NULL, NULL, &post_its));
                PetscCall(PicurvAssertBool((PetscBool)(pre_smoother != post_smoother), "unequal sweeps should give a separate post-smoother"));
                PetscCall(PicurvAssertIntEqual(1, pre_its, "pre-smoothing should run pre_sweeps iterations"));
                PetscCall(PicurvAssertIntEqual(3, post_its, "post-smoothing should run post_sweeps iterations"));
                PetscCall(PetscObjectTypeCompare((PetscObject)post_smoother, KSPRICHARDSON, &is_richardson));
                PetscCall(PicurvAssertBool(is_richardson, "the post-smoother should keep the configured smoother type"));
            } else {
                /* PCMGGetSmootherUp() would itself split a shared smoother, so only the
                   shared one is inspected here. */
                PetscCall(PicurvAssertIntEqual(2, pre_its, "the shared smoother should run the configured sweeps"));
            }
        }
 
        if (c == 0) {
            PetscCall(UpdatePressure(user));
            PetscCall(ProjectVelocity(user));
            PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCONT));
            PetscCall(ComputeDivergence(user));
            PetscCall(PicurvAssertBool((PetscBool)(simCtx->MaxDiv < 1.0e-9 * before),
                                       "the solve with distinct sweeps should still remove the divergence"));
        }
        PetscCall(PicurvDestroyRuntimeContext(&simCtx));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestPoissonNullSpaceRemovesTheInteriorMean()

static PetscErrorCode TestPoissonNullSpaceRemovesTheInteriorMean ( void  )

static

Tests the null space attached to the operator.

Removing it from a constant must leave zero everywhere, and from a linear field must leave zero mean over the interior with zero dummy layers.

Definition at line 845 of file test_poisson_rhs.c.

{
    SimCtx       *simCtx = NULL;
    UserCtx      *user = NULL;
    char          tmpdir[PETSC_MAX_PATH_LEN];
    MatNullSpace  nullsp = NULL;
    Vec           x = NULL;
    PetscReal  ***xa, norm, sum;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(
        NULL, PETSC_FALSE, "-im 9\n-jm 9\n-km 9\n-mg_level 2\n", &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(PoissonSolver_Multigrid(&simCtx->usermg));
    PetscCall(MatGetNullSpace(user->A, &nullsp));
    PetscCall(PicurvAssertBool((PetscBool)(nullsp != NULL), "the operator should carry the Neumann null space"));
 
    PetscCall(VecDuplicate(user->Phi, &x));
    PetscCall(VecSet(x, 3.0));
    PetscCall(MatNullSpaceRemove(nullsp, x));
    PetscCall(VecNorm(x, NORM_INFINITY, &norm));
    PetscCall(PicurvAssertRealNear(0.0, norm, 1.0e-12, "removing the null space should annihilate a constant"));
 
    PetscCall(DMDAVecGetArray(user->da, x, &xa));
    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++) xa[k][j][i] = (PetscReal)i;
    PetscCall(DMDAVecRestoreArray(user->da, x, &xa));
    PetscCall(MatNullSpaceRemove(nullsp, x));
    PetscCall(VecSum(x, &sum));
    PetscCall(PicurvAssertRealNear(0.0, sum, 1.0e-10, "the interior mean should vanish"));
    PetscCall(DMDAVecGetArray(user->da, x, &xa));
    if (user->info.xs == 0) PetscCall(PicurvAssertRealNear(0.0, xa[user->info.zs][user->info.ys][0], 1.0e-14, "dummy layers should be zero"));
    PetscCall(DMDAVecRestoreArray(user->da, x, &xa));
 
    PetscCall(VecDestroy(&x));
    PetscCall(PicurvDestroyRuntimeContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestPoissonSolverMultigridRefusesAnOvercoarsenedHierarchy()

static PetscErrorCode TestPoissonSolverMultigridRefusesAnOvercoarsenedHierarchy ( void  )

static

Tests that an over-coarsened hierarchy stops the run at the Poisson solve.

Three levels on a 9-node grid leave a coarse operator the preconditioner cannot factor. The solve must return PETSC_ERR_NOT_CONVERGED rather than hand an unsolved Phi to the projection.

Definition at line 891 of file test_poisson_rhs.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
    PetscErrorCode solve_error;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(
        NULL, PETSC_FALSE,
        "-im 9\n-jm 9\n-km 9\n-mg_level 3\n"
        /* The shipped profiles solve the coarse level with a replicated direct factor. */
        "-ps_ksp_type fgmres\n-ps_pc_type mg\n"
        "-ps_mg_coarse_ksp_type preonly\n-ps_mg_coarse_pc_type redundant\n",
        &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(PerturbInteriorFaceFluxes(user, 0.2));
 
    PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
    solve_error = PoissonSolver_Multigrid(&simCtx->usermg);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_NOT_CONVERGED, (PetscInt)solve_error,
                                   "an unfactorable coarse level should stop the solve, not reach the projection"));
 
    PetscCall(PicurvDestroyRuntimeContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    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-poisson-rhs PETSc test binary.

Definition at line 922 of file test_poisson_rhs.c.

{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"update-pressure-adds-phi", TestUpdatePressureAddsPhi},
        {"poisson-rhs-zero-divergence", TestPoissonRHSZeroDivergence},
        {"compute-body-forces-dispatcher", TestComputeBodyForcesDispatcher},
        {"compute-eulerian-diffusivity-molecular-only", TestComputeEulerianDiffusivityMolecularOnly},
        {"convection-zero-field", TestConvectionZeroField},
        {"viscous-uniform-field", TestViscousUniformField},
        {"clark-gradient-model-scales-with-filter-width-squared", TestClarkGradientModelScalesWithFilterWidthSquared},
        {"compute-rhs-zero-field-no-forcing", TestComputeRHSZeroFieldNoForcing},
        {"compute-eulerian-diffusivity-gradient-constant-field", TestComputeEulerianDiffusivityGradientConstantField},
        {"compute-eulerian-diffusivity-verification-linear-x", TestComputeEulerianDiffusivityVerificationLinearX},
        {"assemble-poisson-operator-populates-rows", TestAssemblePoissonOperatorPopulatesRows},
        {"operator-and-projection-share-one-face-gradient", TestOperatorAndProjectionShareOneFaceGradient},
        {"projection-zero-phi-leaves-velocity-unchanged", TestProjectionZeroPhiLeavesVelocityUnchanged},
        {"projection-linear-phi-corrects-velocity", TestProjectionLinearPhiCorrectsVelocity},
        {"poisson-solver-multigrid-projects-to-divergence-free", TestPoissonSolverMultigridProjectsToDivergenceFree},
        {"poisson-solver-multigrid-reuses-its-solver", TestPoissonSolverMultigridReusesItsSolver},
        {"poisson-solver-multigrid-honours-distinct-sweeps", TestPoissonSolverMultigridHonoursDistinctSweeps},
        {"poisson-null-space-removes-the-interior-mean", TestPoissonNullSpaceRemovesTheInteriorMean},
        {"poisson-solver-multigrid-refuses-an-overcoarsened-hierarchy", TestPoissonSolverMultigridRefusesAnOvercoarsenedHierarchy},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv Poisson/RHS tests");
    if (ierr) {
        return (int)ierr;
    }
 
    ierr = PicurvRunTests("unit-poisson-rhs", 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: