test_boundaries.c

Go to the documentation of this file.

/**
 * @file test_boundaries.c
 * @brief C unit tests for boundary factories, inlet ownership, and face helpers.
 */
 
#include "test_support.h"
 
#include "Boundaries.h"
#include "wallfunction.h"
 
/**
 * @brief Appends one key/value pair to a linked list of boundary-condition parameters.
 */
static PetscErrorCode AppendBCParam(BC_Param **head, const char *key, const char *value)
{
    BC_Param *node = NULL;
    BC_Param **cursor = head;
 
    PetscFunctionBeginUser;
    PetscCheck(head != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BC param list head cannot be NULL.");
    PetscCall(PetscCalloc1(1, &node));
    PetscCall(PetscStrallocpy(key, &node->key));
    PetscCall(PetscStrallocpy(value, &node->value));
    while (*cursor) {
        cursor = &(*cursor)->next;
    }
    *cursor = node;
    PetscFunctionReturn(0);
}
/**
 * @brief Destroys one boundary-condition handler allocated by a boundary test.
 */
 
static PetscErrorCode DestroyBoundaryHandler(BoundaryCondition **bc_ptr)
{
    BoundaryCondition *bc = NULL;
    PetscFunctionBeginUser;
    if (!bc_ptr || !*bc_ptr) PetscFunctionReturn(0);
 
    bc = *bc_ptr;
    if (bc->Destroy) {
        PetscCall(bc->Destroy(bc));
    }
    PetscCall(PetscFree(bc));
    *bc_ptr = NULL;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Chooses a stable interior node index for face-sample assertions on tiny test grids.
 */
static PetscInt InteriorSampleIndex(PetscInt nodes)
{
    PetscInt index = 3;
 
    if (index >= nodes - 1) {
        index = nodes - 2;
    }
    if (index < 1) {
        index = 1;
    }
    return index;
}
 
/**
 * @brief Extracts the velocity component aligned with the supplied boundary-face normal.
 */
static PetscReal GetFaceNormalComponent(Cmpnts value, BCFace face)
{
    switch (face) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            return value.x;
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            return value.y;
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
            return value.z;
    }
    return 0.0;
}
 
/**
 * @brief Returns the sign convention used for face-normal flux expectations.
 */
static PetscReal GetFaceOrientationSign(BCFace face)
{
    switch (face) {
        case BC_FACE_NEG_X:
        case BC_FACE_NEG_Y:
        case BC_FACE_NEG_Z:
            return 1.0;
        case BC_FACE_POS_X:
        case BC_FACE_POS_Y:
        case BC_FACE_POS_Z:
            return -1.0;
    }
    return 0.0;
}
 
/**
 * @brief Maps an inlet face to the matching configuration key used by the handler parser.
 */
static const char *GetInletParamKey(BCFace face)
{
    switch (face) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            return "vx";
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            return "vy";
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
            return "vz";
    }
    return "";
}
 
/**
 * @brief Computes the number of face-interior sample points for a given boundary face.
 */
static PetscReal GetFaceInteriorPointCount(const UserCtx *user, BCFace face)
{
    switch (face) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            return (PetscReal)(user->info.my - 2) * (PetscReal)(user->info.mz - 2);
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            return (PetscReal)(user->info.mx - 2) * (PetscReal)(user->info.mz - 2);
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
            return (PetscReal)(user->info.mx - 2) * (PetscReal)(user->info.my - 2);
    }
    return 0.0;
}
 
/**
 * @brief Selects representative `Ucont` and `Ubcs` slots for face-matrix assertions.
 */
static PetscErrorCode GetRepresentativeFaceSlots(const UserCtx *user,
                                                 BCFace face,
                                                 PetscInt *ucont_k,
                                                 PetscInt *ucont_j,
                                                 PetscInt *ucont_i,
                                                 PetscInt *ubcs_k,
                                                 PetscInt *ubcs_j,
                                                 PetscInt *ubcs_i)
{
    const PetscInt sample_i = InteriorSampleIndex(user->info.mx);
    const PetscInt sample_j = InteriorSampleIndex(user->info.my);
    const PetscInt sample_k = InteriorSampleIndex(user->info.mz);
 
    PetscFunctionBeginUser;
    switch (face) {
        case BC_FACE_NEG_X:
            *ucont_k = sample_k; *ucont_j = sample_j; *ucont_i = 0;
            *ubcs_k = sample_k; *ubcs_j = sample_j; *ubcs_i = 0;
            break;
        case BC_FACE_POS_X:
            *ucont_k = sample_k; *ucont_j = sample_j; *ucont_i = user->info.mx - 2;
            *ubcs_k = sample_k; *ubcs_j = sample_j; *ubcs_i = user->info.mx - 1;
            break;
        case BC_FACE_NEG_Y:
            *ucont_k = sample_k; *ucont_j = 0; *ucont_i = sample_i;
            *ubcs_k = sample_k; *ubcs_j = 0; *ubcs_i = sample_i;
            break;
        case BC_FACE_POS_Y:
            *ucont_k = sample_k; *ucont_j = user->info.my - 2; *ucont_i = sample_i;
            *ubcs_k = sample_k; *ubcs_j = user->info.my - 1; *ubcs_i = sample_i;
            break;
        case BC_FACE_NEG_Z:
            *ucont_k = 0; *ucont_j = sample_j; *ucont_i = sample_i;
            *ubcs_k = 0; *ubcs_j = sample_j; *ubcs_i = sample_i;
            break;
        case BC_FACE_POS_Z:
            *ucont_k = user->info.mz - 2; *ucont_j = sample_j; *ucont_i = sample_i;
            *ubcs_k = user->info.mz - 1; *ubcs_j = sample_j; *ubcs_i = sample_i;
            break;
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Selects center, off-center, and wall-adjacent slots for parabolic-face checks.
 */
static PetscErrorCode GetParabolicSampleSlots(const UserCtx *user,
                                              BCFace face,
                                              PetscInt *center_k,
                                              PetscInt *center_j,
                                              PetscInt *center_i,
                                              PetscInt *off_k,
                                              PetscInt *off_j,
                                              PetscInt *off_i,
                                              PetscInt *wall_k,
                                              PetscInt *wall_j,
                                              PetscInt *wall_i,
                                              PetscInt *ubcs_center_k,
                                              PetscInt *ubcs_center_j,
                                              PetscInt *ubcs_center_i)
{
    PetscFunctionBeginUser;
    switch (face) {
        case BC_FACE_NEG_X:
            *center_k = 3; *center_j = 3; *center_i = 0;
            *off_k = 2; *off_j = 2; *off_i = 0;
            *wall_k = 1; *wall_j = 1; *wall_i = 0;
            *ubcs_center_k = 3; *ubcs_center_j = 3; *ubcs_center_i = 0;
            break;
        case BC_FACE_POS_X:
            *center_k = 3; *center_j = 3; *center_i = user->info.mx - 2;
            *off_k = 2; *off_j = 2; *off_i = user->info.mx - 2;
            *wall_k = 1; *wall_j = 1; *wall_i = user->info.mx - 2;
            *ubcs_center_k = 3; *ubcs_center_j = 3; *ubcs_center_i = user->info.mx - 1;
            break;
        case BC_FACE_NEG_Y:
            *center_k = 3; *center_j = 0; *center_i = 3;
            *off_k = 2; *off_j = 0; *off_i = 2;
            *wall_k = 1; *wall_j = 0; *wall_i = 1;
            *ubcs_center_k = 3; *ubcs_center_j = 0; *ubcs_center_i = 3;
            break;
        case BC_FACE_POS_Y:
            *center_k = 3; *center_j = user->info.my - 2; *center_i = 3;
            *off_k = 2; *off_j = user->info.my - 2; *off_i = 2;
            *wall_k = 1; *wall_j = user->info.my - 2; *wall_i = 1;
            *ubcs_center_k = 3; *ubcs_center_j = user->info.my - 1; *ubcs_center_i = 3;
            break;
        case BC_FACE_NEG_Z:
            *center_k = 0; *center_j = 3; *center_i = 3;
            *off_k = 0; *off_j = 2; *off_i = 2;
            *wall_k = 0; *wall_j = 1; *wall_i = 1;
            *ubcs_center_k = 0; *ubcs_center_j = 3; *ubcs_center_i = 3;
            break;
        case BC_FACE_POS_Z:
            *center_k = user->info.mz - 2; *center_j = 3; *center_i = 3;
            *off_k = user->info.mz - 2; *off_j = 2; *off_i = 2;
            *wall_k = user->info.mz - 2; *wall_j = 1; *wall_i = 1;
            *ubcs_center_k = user->info.mz - 1; *ubcs_center_j = 3; *ubcs_center_i = 3;
            break;
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Resets one boundary-face configuration entry to a neutral test-local baseline.
 */
static void ResetBoundaryFaceConfig(BoundaryFaceConfig *cfg)
{
    cfg->mathematical_type = INTERFACE;
    cfg->handler_type = BC_HANDLER_UNDEFINED;
    cfg->face_id = BC_FACE_NEG_X;
    if (cfg->params) {
        FreeBC_ParamList(cfg->params);
        cfg->params = NULL;
    }
}
/**
 * @brief Tests that face-service detection matches a defined inlet face.
 */
 
static PetscErrorCode TestCanRankServiceFaceMatchesInletWhenDefined(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscBool generic_face = PETSC_FALSE;
    PetscBool inlet_face = PETSC_FALSE;
    PetscBool any_face_serviceable = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
 
    for (BCFace face = BC_FACE_NEG_X; face <= BC_FACE_POS_Z; face++) {
        PetscCall(CanRankServiceFace(&user->info, user->IM, user->JM, user->KM, face, &generic_face));
        if (generic_face) {
            any_face_serviceable = PETSC_TRUE;
        }
 
        user->inletFaceDefined = PETSC_TRUE;
        user->identifiedInletBCFace = face;
        PetscCall(CanRankServiceInletFace(user, &user->info, user->IM, user->JM, user->KM, &inlet_face));
        PetscCall(PicurvAssertBool((PetscBool)(inlet_face == generic_face), "inlet face service check should match generic face check"));
    }
    PetscCall(PicurvAssertBool(any_face_serviceable, "at least one global face should be serviceable in a non-degenerate domain"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests that inlet-face service requires an inlet face to be defined.
 */
 
static PetscErrorCode TestCanRankServiceInletFaceRequiresDefinition(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscBool can_service = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
 
    user->inletFaceDefined = PETSC_FALSE;
    for (BCFace face = BC_FACE_NEG_X; face <= BC_FACE_POS_Z; face++) {
        user->identifiedInletBCFace = face;
        PetscCall(CanRankServiceInletFace(user, &user->info, user->IM, user->JM, user->KM, &can_service));
        PetscCall(PicurvAssertBool((PetscBool)!can_service, "undefined inlet face should not be serviceable"));
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests the boundary-condition factory assignments for representative handlers.
 */
 
static PetscErrorCode TestBoundaryConditionFactoryAssignments(void)
{
    BoundaryCondition *wall = NULL;
    BoundaryCondition *inlet = NULL;
    BoundaryCondition *parabolic = NULL;
    BoundaryCondition *profile = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_WALL_NOSLIP, &wall));
    PetscCall(PicurvAssertBool((PetscBool)(wall != NULL), "factory should return wall handler"));
    PetscCall(PicurvAssertIntEqual(BC_PRIORITY_WALL, wall->priority, "wall handler priority"));
    PetscCall(PicurvAssertBool((PetscBool)(wall->Apply != NULL), "wall handler should expose Apply"));
    PetscCall(PicurvAssertBool((PetscBool)(wall->Destroy == NULL), "wall handler should not require destroy hook"));
    PetscCall(DestroyBoundaryHandler(&wall));
 
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_CONSTANT_VELOCITY, &inlet));
    PetscCall(PicurvAssertBool((PetscBool)(inlet != NULL), "factory should return inlet handler"));
    PetscCall(PicurvAssertIntEqual(BC_PRIORITY_INLET, inlet->priority, "inlet handler priority"));
    PetscCall(PicurvAssertBool((PetscBool)(inlet->Initialize != NULL), "inlet handler should expose Initialize"));
    PetscCall(PicurvAssertBool((PetscBool)(inlet->Apply != NULL), "inlet handler should expose Apply"));
    PetscCall(PicurvAssertBool((PetscBool)(inlet->Destroy != NULL), "inlet handler should expose its destroy hook"));
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PARABOLIC, &parabolic));
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PROFILE_FROM_FILE, &profile));
    PetscCall(PicurvAssertBool((PetscBool)(inlet->Apply == parabolic->Apply &&
                                           inlet->Apply == profile->Apply),
                               "all profile-producing inlet modes should share one application hook"));
    PetscCall(DestroyBoundaryHandler(&inlet));
    PetscCall(DestroyBoundaryHandler(&parabolic));
    PetscCall(DestroyBoundaryHandler(&profile));
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Reapplies a static inlet through the first three time-step contexts.
 */
static PetscErrorCode ApplyInletAcrossFirstSteps(BoundaryCondition *bc, BCContext *ctx)
{
    PetscFunctionBeginUser;
    for (PetscInt step = 0; step < 3; step++) {
        PetscReal local_inflow = 0.0, local_outflow = 0.0;
        ctx->user->simCtx->step = step;
        ctx->user->simCtx->ti = 0.1 * (PetscReal)step;
        PetscCall(VecSet(ctx->user->Ucont, -91.0));
        PetscCall(VecSet(ctx->user->Bcs.Ubcs, -92.0));
        PetscCall(bc->PreStep(bc, ctx, &local_inflow, &local_outflow));
        PetscCall(bc->Apply(bc, ctx));
    }
    PetscFunctionReturn(0);
}
/**
 * @brief Tests the implemented-handler matrix exposed by the factory.
 */
 
static PetscErrorCode TestBoundaryConditionFactoryImplementedHandlerMatrix(void)
{
    struct HandlerExpectation {
        BCHandlerType handler;
        BCPriorityType priority;
        PetscBool expect_apply;
        PetscBool expect_initialize;
    };
    const struct HandlerExpectation expectations[] = {
        {BC_HANDLER_WALL_NOSLIP, BC_PRIORITY_WALL, PETSC_TRUE, PETSC_FALSE},
        {BC_HANDLER_INLET_CONSTANT_VELOCITY, BC_PRIORITY_INLET, PETSC_TRUE, PETSC_TRUE},
        {BC_HANDLER_INLET_PARABOLIC, BC_PRIORITY_INLET, PETSC_TRUE, PETSC_TRUE},
        {BC_HANDLER_INLET_PROFILE_FROM_FILE, BC_PRIORITY_INLET, PETSC_TRUE, PETSC_TRUE},
        {BC_HANDLER_OUTLET_CONSERVATION, BC_PRIORITY_OUTLET, PETSC_TRUE, PETSC_FALSE},
    };
 
    PetscFunctionBeginUser;
    for (size_t i = 0; i < sizeof(expectations) / sizeof(expectations[0]); ++i) {
        BoundaryCondition *bc = NULL;
        PetscCall(BoundaryCondition_Create(expectations[i].handler, &bc));
        PetscCall(PicurvAssertBool((PetscBool)(bc != NULL), "factory should allocate a handler object"));
        PetscCall(PicurvAssertIntEqual(expectations[i].priority, bc->priority, "handler priority should match expectation"));
        PetscCall(PicurvAssertBool((PetscBool)((bc->Apply != NULL) == expectations[i].expect_apply), "Apply hook expectation mismatch"));
        PetscCall(PicurvAssertBool((PetscBool)((bc->Initialize != NULL) == expectations[i].expect_initialize), "Initialize hook expectation mismatch"));
        PetscCall(DestroyBoundaryHandler(&bc));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Writes a tiny canonical PICSLICE profile for boundary handler tests.
 */
static PetscErrorCode WritePicSliceForTests(const char *path, PetscInt n1, PetscInt n2, PetscReal base)
{
    FILE *fd = NULL;
 
    PetscFunctionBeginUser;
    fd = fopen(path, "w");
    PetscCheck(fd != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open %s for writing.", path);
    PetscCheck(fprintf(fd, "PICSLICE\n1\n%d %d\n", (int)n1, (int)n2) >= 0,
               PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE, "Could not write PICSLICE header.");
    for (PetscInt a = 0; a < n1; a++) {
        for (PetscInt b = 0; b < n2; b++) {
            PetscReal value = base + (PetscReal)(10 * a + b);
            PetscCheck(fprintf(fd, "%.16e\n", (double)value) >= 0,
                       PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE, "Could not write PICSLICE value.");
        }
    }
    fclose(fd);
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests prescribed-flow inlet handler loading and applying one Z-face profile.
 */
static PetscErrorCode TestInletProfileFromFileHandlerBehavior(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
    char profile_path[PETSC_MAX_PATH_LEN];
    PetscReal local_inflow = 0.0;
    PetscReal local_outflow = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
    PetscCall(PicurvPopulateIdentityMetrics(user));
 
    PetscCall(PetscSNPrintf(profile_path, sizeof(profile_path), "/tmp/picurv_unit_profile_%d.picslice", (int)getpid()));
    PetscCall(WritePicSliceForTests(profile_path, user->JM - 1, user->IM - 1, 2.0));
 
    user->boundary_faces[BC_FACE_NEG_Z].face_id = BC_FACE_NEG_Z;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = INLET;
    user->boundary_faces[BC_FACE_NEG_Z].handler_type = BC_HANDLER_INLET_PROFILE_FROM_FILE;
    PetscCall(AppendBCParam(&user->boundary_faces[BC_FACE_NEG_Z].params, "source_file", profile_path));
 
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
 
    ctx.user = user;
    ctx.face_id = BC_FACE_NEG_Z;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PROFILE_FROM_FILE, &bc));
    PetscCall(bc->PreStep(bc, &ctx, &local_inflow, &local_outflow));
    PetscCall(PicurvAssertRealNear(0.0, local_inflow, 1.0e-12, "profile inlet PreStep should leave inflow unchanged"));
    PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "profile inlet PreStep should leave outflow unchanged"));
    PetscCall(bc->Initialize(bc, &ctx));
    PetscCall(ApplyInletAcrossFirstSteps(bc, &ctx));
    PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(PicurvAssertRealNear(2.0, ucont[0][1][1].z, 1.0e-12, "profile inlet should apply first scalar speed"));
    PetscCall(PicurvAssertRealNear(13.0, ucont[0][2][2].z, 1.0e-12, "profile inlet should preserve PICSLICE ordering"));
    PetscCall(PicurvAssertRealNear(24.0, ubcs[0][3][3].z, 1.0e-12, "profile inlet should write boundary velocity"));
    PetscCall(PicurvAssertRealNear(35.0, ucont[0][4][4].z, 1.0e-12, "profile inlet should consume the last face-cell slot"));
    PetscCall(PicurvAssertRealNear(35.0, ubcs[0][4][4].z, 1.0e-12, "profile inlet should write the last face-cell boundary value"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(PicurvAssertBool((PetscBool)(local_inflow > 0.0), "profile inlet PostStep should accumulate positive negative-face flux"));
    PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "profile inlet should not add to outflow"));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    FreeBC_ParamList(user->boundary_faces[BC_FACE_NEG_Z].params);
    user->boundary_faces[BC_FACE_NEG_Z].params = NULL;
    remove(profile_path);
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests that unsupported handler types are rejected by the factory.
 */
 
static PetscErrorCode TestBoundaryConditionFactoryRejectsUnsupportedHandler(void)
{
    BoundaryCondition *bc = NULL;
    PetscErrorCode ierr_create;
 
    PetscFunctionBeginUser;
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    ierr_create = BoundaryCondition_Create(BC_HANDLER_OUTLET_PRESSURE, &bc);
    PetscCall(PetscPopErrorHandler());
 
    PetscCall(PicurvAssertBool((PetscBool)(ierr_create != 0), "unsupported handler should return a non-zero error code"));
    if (bc) {
        PetscCall(PetscFree(bc));
    }
    PetscFunctionReturn(0);
}
/**
 * @brief Tests deterministic inlet-face grid-location helpers across all faces.
 */
 
static PetscErrorCode TestGetDeterministicFaceGridLocationFaceMatrix(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 8, 8));
    simCtx->np = 128;
 
    for (BCFace face = BC_FACE_NEG_X; face <= BC_FACE_POS_Z; ++face) {
        PetscInt ci = -1, cj = -1, ck = -1;
        PetscReal xi = -1.0, eta = -1.0, zta = -1.0;
        PetscBool placed = PETSC_FALSE;
 
        user->identifiedInletBCFace = face;
        PetscCall(GetDeterministicFaceGridLocation(
            user, &user->info,
            user->info.xs, user->info.ys, user->info.zs,
            user->info.mx, user->info.my, user->info.mz,
            0,
            &ci, &cj, &ck, &xi, &eta, &zta, &placed));
 
        PetscCall(PicurvAssertBool(placed, "single-rank deterministic face placement should succeed"));
        PetscCall(PicurvAssertBool((PetscBool)(ci >= user->info.xs && ci < user->info.xs + user->info.xm), "ci must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(cj >= user->info.ys && cj < user->info.ys + user->info.ym), "cj must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(ck >= user->info.zs && ck < user->info.zs + user->info.zm), "ck must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(xi >= 0.0 && xi < 1.0), "xi should be in [0,1)"));
        PetscCall(PicurvAssertBool((PetscBool)(eta >= 0.0 && eta < 1.0), "eta should be in [0,1)"));
        PetscCall(PicurvAssertBool((PetscBool)(zta >= 0.0 && zta < 1.0), "zta should be in [0,1)"));
 
        if (face == BC_FACE_NEG_X || face == BC_FACE_POS_X) {
            PetscCall(PicurvAssertRealNear(0.5, xi, 1.0e-10, "deterministic x-face placement should sit halfway into boundary-adjacent cell"));
        } else if (face == BC_FACE_NEG_Y || face == BC_FACE_POS_Y) {
            PetscCall(PicurvAssertRealNear(0.5, eta, 1.0e-10, "deterministic y-face placement should sit halfway into boundary-adjacent cell"));
        } else {
            PetscCall(PicurvAssertRealNear(0.5, zta, 1.0e-10, "deterministic z-face placement should sit halfway into boundary-adjacent cell"));
        }
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests random inlet-face cell selection across all supported faces.
 */
 
static PetscErrorCode TestGetRandomCellAndLogicalCoordsOnInletFaceMatrix(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscRandom rand_i = NULL, rand_j = NULL, rand_k = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 8, 8));
 
    PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand_i));
    PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand_j));
    PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand_k));
    PetscCall(PetscRandomSetInterval(rand_i, 0.0, 1.0));
    PetscCall(PetscRandomSetInterval(rand_j, 0.0, 1.0));
    PetscCall(PetscRandomSetInterval(rand_k, 0.0, 1.0));
    PetscCall(PetscRandomSetFromOptions(rand_i));
    PetscCall(PetscRandomSetFromOptions(rand_j));
    PetscCall(PetscRandomSetFromOptions(rand_k));
 
    for (BCFace face = BC_FACE_NEG_X; face <= BC_FACE_POS_Z; ++face) {
        PetscInt ci = -1, cj = -1, ck = -1;
        PetscReal xi = -1.0, eta = -1.0, zta = -1.0;
        const PetscReal boundary_eps = 5.0e-4;
 
        user->identifiedInletBCFace = face;
        PetscCall(GetRandomCellAndLogicalCoordsOnInletFace(
            user, &user->info,
            user->info.xs, user->info.ys, user->info.zs,
            user->info.mx, user->info.my, user->info.mz,
            &rand_i, &rand_j, &rand_k,
            &ci, &cj, &ck,
            &xi, &eta, &zta));
 
        PetscCall(PicurvAssertBool((PetscBool)(ci >= user->info.xs && ci < user->info.xs + user->info.xm), "ci must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(cj >= user->info.ys && cj < user->info.ys + user->info.ym), "cj must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(ck >= user->info.zs && ck < user->info.zs + user->info.zm), "ck must map to owned node window"));
        PetscCall(PicurvAssertBool((PetscBool)(xi >= 0.0 && xi <= 1.0), "xi should be in [0,1]"));
        PetscCall(PicurvAssertBool((PetscBool)(eta >= 0.0 && eta <= 1.0), "eta should be in [0,1]"));
        PetscCall(PicurvAssertBool((PetscBool)(zta >= 0.0 && zta <= 1.0), "zta should be in [0,1]"));
 
        switch (face) {
            case BC_FACE_NEG_X:
                PetscCall(PicurvAssertBool((PetscBool)(xi <= boundary_eps), "NEG_X should pin xi near 0"));
                break;
            case BC_FACE_POS_X:
                PetscCall(PicurvAssertBool((PetscBool)(xi >= 1.0 - boundary_eps), "POS_X should pin xi near 1"));
                break;
            case BC_FACE_NEG_Y:
                PetscCall(PicurvAssertBool((PetscBool)(eta <= boundary_eps), "NEG_Y should pin eta near 0"));
                break;
            case BC_FACE_POS_Y:
                PetscCall(PicurvAssertBool((PetscBool)(eta >= 1.0 - boundary_eps), "POS_Y should pin eta near 1"));
                break;
            case BC_FACE_NEG_Z:
                PetscCall(PicurvAssertBool((PetscBool)(zta <= boundary_eps), "NEG_Z should pin zta near 0"));
                break;
            case BC_FACE_POS_Z:
                PetscCall(PicurvAssertBool((PetscBool)(zta >= 1.0 - boundary_eps), "POS_Z should pin zta near 1"));
                break;
        }
    }
 
    PetscCall(PetscRandomDestroy(&rand_i));
    PetscCall(PetscRandomDestroy(&rand_j));
    PetscCall(PetscRandomDestroy(&rand_k));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests constant inlet handler initialization and face writes on a Z inlet.
 */
 
static PetscErrorCode TestInletConstantVelocityHandlerBehavior(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
    PetscReal local_inflow = 0.0;
    PetscReal local_outflow = 0.0;
    PetscReal expected_flux = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
    PetscCall(PicurvPopulateIdentityMetrics(user));
 
    user->boundary_faces[BC_FACE_NEG_Z].face_id = BC_FACE_NEG_Z;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = INLET;
    user->boundary_faces[BC_FACE_NEG_Z].handler_type = BC_HANDLER_INLET_CONSTANT_VELOCITY;
    PetscCall(AppendBCParam(&user->boundary_faces[BC_FACE_NEG_Z].params, "vz", "2.5"));
 
    ctx.user = user;
    ctx.face_id = BC_FACE_NEG_Z;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_CONSTANT_VELOCITY, &bc));
    PetscCall(bc->Initialize(bc, &ctx));
    PetscCall(ApplyInletAcrossFirstSteps(bc, &ctx));
    PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(PicurvAssertRealNear(2.5, ucont[0][3][3].z, 1.0e-12, "constant inlet should set Ucont normal component"));
    PetscCall(PicurvAssertRealNear(2.5, ubcs[0][3][3].z, 1.0e-12, "constant inlet should set Ubcs normal component"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
 
    expected_flux = (PetscReal)(user->info.mx - 2) * (PetscReal)(user->info.my - 2) * 2.5;
    PetscCall(PicurvAssertRealNear(expected_flux, local_inflow, 1.0e-12, "constant inlet PostStep should sum the face flux"));
    PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "constant inlet should not contribute to outflow accumulation"));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    FreeBC_ParamList(user->boundary_faces[BC_FACE_NEG_Z].params);
    user->boundary_faces[BC_FACE_NEG_Z].params = NULL;
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Verifies curvilinear projection and Cartesian dummy reflection for the shared inlet path.
 */
static PetscErrorCode TestInletCommonApplicationCurvilinearMetric(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***zet = NULL, ***ucont = NULL, ***ubcs = NULL, ***ucat = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
    PetscCall(DMDAVecGetArray(user->fda, user->lZet, &zet));
    zet[0][3][3].x = 0.0;
    zet[0][3][3].y = 3.0;
    zet[0][3][3].z = 4.0;
    PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &zet));
 
    user->boundary_faces[BC_FACE_NEG_Z].face_id = BC_FACE_NEG_Z;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = INLET;
    user->boundary_faces[BC_FACE_NEG_Z].handler_type = BC_HANDLER_INLET_CONSTANT_VELOCITY;
    PetscCall(AppendBCParam(&user->boundary_faces[BC_FACE_NEG_Z].params, "vz", "2.0"));
    ctx.user = user;
    ctx.face_id = BC_FACE_NEG_Z;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_CONSTANT_VELOCITY, &bc));
    PetscCall(bc->Initialize(bc, &ctx));
    PetscCall(bc->Apply(bc, &ctx));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    PetscCall(PicurvAssertRealNear(10.0, ucont[0][3][3].z, 1.0e-12,
                                   "Ucont should be the Cartesian inlet velocity dotted with the face metric"));
    PetscCall(PicurvAssertRealNear(0.0, ubcs[0][3][3].x, 1.0e-12,
                                   "curvilinear normal velocity should preserve its Cartesian x component"));
    PetscCall(PicurvAssertRealNear(1.2, ubcs[0][3][3].y, 1.0e-12,
                                   "curvilinear normal velocity should include Cartesian y"));
    PetscCall(PicurvAssertRealNear(1.6, ubcs[0][3][3].z, 1.0e-12,
                                   "curvilinear normal velocity should include Cartesian z"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
 
    PetscCall(VecSet(user->Ucat, 1.0));
    PetscCall(UpdateDummyCells(user, FIELD_ID_UCAT));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(PicurvAssertRealNear(-1.0, ucat[0][3][3].x, 1.0e-12,
                                   "dummy reflection should act on Cartesian x"));
    PetscCall(PicurvAssertRealNear(1.4, ucat[0][3][3].y, 1.0e-12,
                                   "dummy reflection should act on Cartesian y"));
    PetscCall(PicurvAssertRealNear(2.2, ucat[0][3][3].z, 1.0e-12,
                                   "dummy reflection should act on Cartesian z"));
    PetscCall(PicurvAssertRealNear(1.0, ucat[1][3][3].y, 1.0e-12,
                                   "dummy reflection should not overwrite the adjacent interior value"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    ResetBoundaryFaceConfig(&user->boundary_faces[BC_FACE_NEG_Z]);
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests wall no-slip application across the full face matrix.
 */
static PetscErrorCode TestWallNoSlipHandlerFaceMatrix(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
    const BCFace faces[] = {
        BC_FACE_NEG_X, BC_FACE_POS_X,
        BC_FACE_NEG_Y, BC_FACE_POS_Y,
        BC_FACE_NEG_Z, BC_FACE_POS_Z
    };
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(BoundaryCondition_Create(BC_HANDLER_WALL_NOSLIP, &bc));
 
    for (size_t n = 0; n < sizeof(faces) / sizeof(faces[0]); ++n) {
        const BCFace face = faces[n];
        PetscInt ucont_k, ucont_j, ucont_i, ubcs_k, ubcs_j, ubcs_i;
 
        PetscCall(VecSet(user->Ucont, 5.0));
        PetscCall(VecSet(user->Bcs.Ubcs, 7.0));
        ctx.user = user;
        ctx.face_id = face;
        PetscCall(bc->Apply(bc, &ctx));
 
        PetscCall(GetRepresentativeFaceSlots(user, face, &ucont_k, &ucont_j, &ucont_i, &ubcs_k, &ubcs_j, &ubcs_i));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(0.0, GetFaceNormalComponent(ucont[ucont_k][ucont_j][ucont_i], face), 1.0e-12, "wall no-slip should zero the normal flux component"));
        PetscCall(PicurvAssertRealNear(0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].x, 1.0e-12, "wall no-slip should zero Ubcs.x"));
        PetscCall(PicurvAssertRealNear(0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].y, 1.0e-12, "wall no-slip should zero Ubcs.y"));
        PetscCall(PicurvAssertRealNear(0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].z, 1.0e-12, "wall no-slip should zero Ubcs.z"));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
    }
 
    PetscCall(DestroyBoundaryHandler(&bc));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests constant-inlet initialization, flux accounting, and face writes across all faces.
 */
static PetscErrorCode TestInletConstantVelocityHandlerFaceMatrix(void)
{
    struct ConstantInletCase {
        BCFace face;
        PetscReal velocity;
        const char *value_text;
    };
    const struct ConstantInletCase cases[] = {
        {BC_FACE_NEG_X, 1.5, "1.5"},
        {BC_FACE_POS_X, 2.0, "2.0"},
        {BC_FACE_NEG_Y, 2.5, "2.5"},
        {BC_FACE_POS_Y, 3.0, "3.0"},
        {BC_FACE_NEG_Z, 3.5, "3.5"},
        {BC_FACE_POS_Z, 4.0, "4.0"},
    };
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
 
    for (size_t n = 0; n < sizeof(cases) / sizeof(cases[0]); ++n) {
        const struct ConstantInletCase *test_case = &cases[n];
        BoundaryCondition *bc = NULL;
        PetscReal local_inflow = 0.0;
        PetscReal local_outflow = 0.0;
        PetscInt ucont_k, ucont_j, ucont_i, ubcs_k, ubcs_j, ubcs_i;
        const PetscReal expected_value = GetFaceOrientationSign(test_case->face) * test_case->velocity;
        const PetscReal expected_flux = GetFaceInteriorPointCount(user, test_case->face) * expected_value;
 
        PetscCall(VecSet(user->Ucont, 0.0));
        PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
        ResetBoundaryFaceConfig(&user->boundary_faces[test_case->face]);
        user->boundary_faces[test_case->face].face_id = test_case->face;
        user->boundary_faces[test_case->face].mathematical_type = INLET;
        user->boundary_faces[test_case->face].handler_type = BC_HANDLER_INLET_CONSTANT_VELOCITY;
        PetscCall(AppendBCParam(&user->boundary_faces[test_case->face].params,
                                GetInletParamKey(test_case->face),
                                test_case->value_text));
 
        ctx.user = user;
        ctx.face_id = test_case->face;
        PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_CONSTANT_VELOCITY, &bc));
        PetscCall(bc->PreStep(bc, &ctx, &local_inflow, &local_outflow));
        PetscCall(PicurvAssertRealNear(0.0, local_inflow, 1.0e-12, "constant inlet PreStep should leave inflow unchanged"));
        PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "constant inlet PreStep should leave outflow unchanged"));
        PetscCall(bc->Initialize(bc, &ctx));
        PetscCall(bc->Apply(bc, &ctx));
        PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
        PetscCall(GetRepresentativeFaceSlots(user, test_case->face, &ucont_k, &ucont_j, &ucont_i, &ubcs_k, &ubcs_j, &ubcs_i));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(expected_value, GetFaceNormalComponent(ucont[ucont_k][ucont_j][ucont_i], test_case->face), 1.0e-12, "constant inlet should set the face-normal Ucont component"));
        PetscCall(PicurvAssertRealNear(expected_value, GetFaceNormalComponent(ubcs[ubcs_k][ubcs_j][ubcs_i], test_case->face), 1.0e-12, "constant inlet should set the face-normal Ubcs component"));
        PetscCall(PicurvAssertRealNear((test_case->face == BC_FACE_NEG_X || test_case->face == BC_FACE_POS_X) ? expected_value : 0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].x, 1.0e-12, "constant inlet should only set the expected Ubcs axis"));
        PetscCall(PicurvAssertRealNear((test_case->face == BC_FACE_NEG_Y || test_case->face == BC_FACE_POS_Y) ? expected_value : 0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].y, 1.0e-12, "constant inlet should only set the expected Ubcs axis"));
        PetscCall(PicurvAssertRealNear((test_case->face == BC_FACE_NEG_Z || test_case->face == BC_FACE_POS_Z) ? expected_value : 0.0, ubcs[ubcs_k][ubcs_j][ubcs_i].z, 1.0e-12, "constant inlet should only set the expected Ubcs axis"));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(expected_flux, local_inflow, 1.0e-12, "constant inlet PostStep should integrate the face flux with orientation"));
        PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "constant inlet should not add to outflow"));
 
        PetscCall(DestroyBoundaryHandler(&bc));
        ResetBoundaryFaceConfig(&user->boundary_faces[test_case->face]);
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests parabolic inlet handler shape on a tiny Z-face.
 */
 
static PetscErrorCode TestInletParabolicProfileHandlerBehavior(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
 
    user->boundary_faces[BC_FACE_NEG_Z].face_id = BC_FACE_NEG_Z;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = INLET;
    user->boundary_faces[BC_FACE_NEG_Z].handler_type = BC_HANDLER_INLET_PARABOLIC;
    PetscCall(AppendBCParam(&user->boundary_faces[BC_FACE_NEG_Z].params, "v_max", "4.0"));
 
    ctx.user = user;
    ctx.face_id = BC_FACE_NEG_Z;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PARABOLIC, &bc));
    PetscCall(bc->Initialize(bc, &ctx));
    PetscCall(ApplyInletAcrossFirstSteps(bc, &ctx));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(PicurvAssertRealNear(4.0, ucont[0][3][3].z, 1.0e-12, "parabolic inlet should peak at the face centerline"));
    /* Five cells per axis: the wall-adjacent centre sits at s = -0.8, so each factor is 0.36. */
    PetscCall(PicurvAssertRealNear(4.0 * 0.36 * 0.36, ucont[0][1][1].z, 1.0e-12,
                                   "parabolic inlet should vanish on the wall, not at the wall-adjacent centre"));
    PetscCall(PicurvAssertBool((PetscBool)(ucont[0][3][3].z > ucont[0][2][2].z), "parabolic inlet centerline should exceed off-center velocity"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    FreeBC_ParamList(user->boundary_faces[BC_FACE_NEG_Z].params);
    user->boundary_faces[BC_FACE_NEG_Z].params = NULL;
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the parabolic inlet is uniform along a periodic cross-stream axis.
 * @details A channel periodic in i has walls only in j, so its inlet is the exact
 *          one-dimensional parabola in j.
 */
static PetscErrorCode TestInletParabolicProfileUniformAlongPeriodicAxis(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
    simCtx->i_periodic = 1;
 
    user->boundary_faces[BC_FACE_NEG_Z].face_id = BC_FACE_NEG_Z;
    user->boundary_faces[BC_FACE_NEG_Z].mathematical_type = INLET;
    user->boundary_faces[BC_FACE_NEG_Z].handler_type = BC_HANDLER_INLET_PARABOLIC;
    PetscCall(AppendBCParam(&user->boundary_faces[BC_FACE_NEG_Z].params, "v_max", "4.0"));
    ctx.user = user;
    ctx.face_id = BC_FACE_NEG_Z;
    PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PARABOLIC, &bc));
    PetscCall(bc->Initialize(bc, &ctx));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(PicurvAssertRealNear(4.0, ucont[0][3][1].z, 1.0e-12, "periodic axis carries no parabola factor"));
    PetscCall(PicurvAssertRealNear(ucont[0][3][1].z, ucont[0][3][4].z, 1.0e-12, "profile is uniform along the periodic axis"));
    PetscCall(PicurvAssertRealNear(4.0 * 0.36, ucont[0][1][3].z, 1.0e-12, "wall-bounded axis keeps its parabola"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    FreeBC_ParamList(user->boundary_faces[BC_FACE_NEG_Z].params);
    user->boundary_faces[BC_FACE_NEG_Z].params = NULL;
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests parabolic-inlet centerline, wall, and flux behavior across all faces.
 */
static PetscErrorCode TestInletParabolicProfileHandlerFaceMatrix(void)
{
    const BCFace faces[] = {
        BC_FACE_NEG_X, BC_FACE_POS_X,
        BC_FACE_NEG_Y, BC_FACE_POS_Y,
        BC_FACE_NEG_Z, BC_FACE_POS_Z
    };
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
 
    for (size_t n = 0; n < sizeof(faces) / sizeof(faces[0]); ++n) {
        const BCFace face = faces[n];
        BoundaryCondition *bc = NULL;
        PetscReal local_inflow = 0.0;
        PetscReal local_outflow = 0.0;
        PetscInt center_k, center_j, center_i;
        PetscInt off_k, off_j, off_i;
        PetscInt wall_k, wall_j, wall_i;
        PetscInt ubcs_center_k, ubcs_center_j, ubcs_center_i;
        const PetscReal expected_center = GetFaceOrientationSign(face) * 4.0;
 
        PetscCall(VecSet(user->Ucont, 0.0));
        PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
        ResetBoundaryFaceConfig(&user->boundary_faces[face]);
        user->boundary_faces[face].face_id = face;
        user->boundary_faces[face].mathematical_type = INLET;
        user->boundary_faces[face].handler_type = BC_HANDLER_INLET_PARABOLIC;
        PetscCall(AppendBCParam(&user->boundary_faces[face].params, "v_max", "4.0"));
 
        ctx.user = user;
        ctx.face_id = face;
        PetscCall(BoundaryCondition_Create(BC_HANDLER_INLET_PARABOLIC, &bc));
        PetscCall(bc->PreStep(bc, &ctx, &local_inflow, &local_outflow));
        PetscCall(PicurvAssertRealNear(0.0, local_inflow, 1.0e-12, "parabolic inlet PreStep should leave inflow unchanged"));
        PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "parabolic inlet PreStep should leave outflow unchanged"));
        PetscCall(bc->Initialize(bc, &ctx));
        PetscCall(bc->Apply(bc, &ctx));
        PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
        PetscCall(GetParabolicSampleSlots(user, face,
                                          &center_k, &center_j, &center_i,
                                          &off_k, &off_j, &off_i,
                                          &wall_k, &wall_j, &wall_i,
                                          &ubcs_center_k, &ubcs_center_j, &ubcs_center_i));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(expected_center, GetFaceNormalComponent(ucont[center_k][center_j][center_i], face), 1.0e-12, "parabolic inlet should peak at the face centerline"));
        PetscCall(PicurvAssertRealNear(expected_center, GetFaceNormalComponent(ubcs[ubcs_center_k][ubcs_center_j][ubcs_center_i], face), 1.0e-12, "parabolic inlet should write the centerline boundary velocity"));
        PetscCall(PicurvAssertRealNear(expected_center * 0.36 * 0.36, GetFaceNormalComponent(ucont[wall_k][wall_j][wall_i], face), 1.0e-12,
                                       "parabolic inlet should take its wall-adjacent value, vanishing only on the wall"));
        PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(GetFaceNormalComponent(ucont[center_k][center_j][center_i], face)) >
                                               PetscAbsReal(GetFaceNormalComponent(ucont[off_k][off_j][off_i], face))),
                                   "parabolic inlet centerline magnitude should exceed the off-center magnitude"));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(local_inflow) > 0.0), "parabolic inlet PostStep should accumulate a non-zero flux"));
        PetscCall(PicurvAssertBool((PetscBool)((local_inflow > 0.0) == (GetFaceOrientationSign(face) > 0.0)),
                                   "parabolic inlet PostStep should preserve the face orientation"));
        PetscCall(PicurvAssertRealNear(0.0, local_outflow, 1.0e-12, "parabolic inlet should not add to outflow"));
 
        PetscCall(DestroyBoundaryHandler(&bc));
        ResetBoundaryFaceConfig(&user->boundary_faces[face]);
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief Tests outlet conservation handler correction and post-step flux accounting.
 */
 
static PetscErrorCode TestOutletConservationHandlerBehavior(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BoundaryCondition *bc = NULL;
    BCContext ctx;
    PetscReal global_inflow = 30.0;
    PetscReal global_farfield_in = 0.0;
    PetscReal global_farfield_out = 0.0;
    PetscReal global_outflow = 25.0;
    PetscReal local_inflow = 0.0;
    PetscReal local_outflow = 0.0;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ucat = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
    PetscCall(VecSet(user->Ucat, 1.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    simCtx->AreaOutSum = (PetscReal)(user->info.mx - 2) * (PetscReal)(user->info.my - 2);
 
    user->boundary_faces[BC_FACE_POS_Z].face_id = BC_FACE_POS_Z;
    user->boundary_faces[BC_FACE_POS_Z].mathematical_type = OUTLET;
    user->boundary_faces[BC_FACE_POS_Z].handler_type = BC_HANDLER_OUTLET_CONSERVATION;
 
    ctx.user = user;
    ctx.face_id = BC_FACE_POS_Z;
    ctx.global_inflow_sum = &global_inflow;
    ctx.global_farfield_inflow_sum = &global_farfield_in;
    ctx.global_farfield_outflow_sum = &global_farfield_out;
    ctx.global_outflow_sum = &global_outflow;
 
    PetscCall(BoundaryCondition_Create(BC_HANDLER_OUTLET_CONSERVATION, &bc));
    PetscCall(bc->Apply(bc, &ctx));
    PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lUcat, &ucat));
    PetscCall(PicurvAssertRealNear(1.2, ucont[5][3][3].z, 1.0e-12, "outlet correction should add the expected flux trim"));
    PetscCall(PicurvAssertRealNear(1.0, ucat[5][3][3].z, 1.0e-12, "outlet handler should preserve the interior Ucat reference"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lUcat, &ucat));
 
    PetscCall(PicurvAssertRealNear(30.0, local_outflow, 1.0e-12, "outlet PostStep should report corrected flux"));
    PetscCall(PicurvAssertRealNear(0.0, local_inflow, 1.0e-12, "outlet handler should not accumulate inflow"));
 
    PetscCall(DestroyBoundaryHandler(&bc));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests outlet-conservation correction and flux accounting across all outlet faces.
 */
static PetscErrorCode TestOutletConservationHandlerFaceMatrix(void)
{
    const BCFace faces[] = {
        BC_FACE_NEG_X, BC_FACE_POS_X,
        BC_FACE_NEG_Y, BC_FACE_POS_Y,
        BC_FACE_NEG_Z, BC_FACE_POS_Z
    };
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    BCContext ctx;
    Cmpnts ***ucont = NULL;
    Cmpnts ***ubcs = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(PicurvPopulateIdentityMetrics(user));
    PetscCall(VecSet(user->Ucat, 1.0));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
 
    for (size_t n = 0; n < sizeof(faces) / sizeof(faces[0]); ++n) {
        const BCFace face = faces[n];
        BoundaryCondition *bc = NULL;
        PetscReal global_inflow = 0.0;
        PetscReal global_farfield_in = 0.0;
        PetscReal global_farfield_out = 0.0;
        PetscReal global_outflow = 0.0;
        PetscReal local_inflow = 0.0;
        PetscReal local_outflow = 0.0;
        PetscReal area = GetFaceInteriorPointCount(user, face);
        PetscInt ucont_k, ucont_j, ucont_i, ubcs_k, ubcs_j, ubcs_i;
 
        PetscCall(VecSet(user->Ucont, 0.0));
        PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
        ResetBoundaryFaceConfig(&user->boundary_faces[face]);
        user->boundary_faces[face].face_id = face;
        user->boundary_faces[face].mathematical_type = OUTLET;
        user->boundary_faces[face].handler_type = BC_HANDLER_OUTLET_CONSERVATION;
        user->simCtx->AreaOutSum = area;
 
        ctx.user = user;
        ctx.face_id = face;
        ctx.global_inflow_sum = &global_inflow;
        ctx.global_farfield_inflow_sum = &global_farfield_in;
        ctx.global_farfield_outflow_sum = &global_farfield_out;
        ctx.global_outflow_sum = &global_outflow;
 
        PetscCall(BoundaryCondition_Create(BC_HANDLER_OUTLET_CONSERVATION, &bc));
        PetscCall(bc->PreStep(bc, &ctx, &local_inflow, &local_outflow));
        PetscCall(PicurvAssertRealNear(area, local_outflow, 1.0e-12, "outlet PreStep should measure the uncorrected face flux"));
 
        global_inflow = area + 0.2 * area;
        global_outflow = local_outflow;
        local_outflow = 0.0;
        PetscCall(bc->Apply(bc, &ctx));
        PetscCall(bc->PostStep(bc, &ctx, &local_inflow, &local_outflow));
 
        PetscCall(GetRepresentativeFaceSlots(user, face, &ucont_k, &ucont_j, &ucont_i, &ubcs_k, &ubcs_j, &ubcs_i));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(1.2, GetFaceNormalComponent(ucont[ucont_k][ucont_j][ucont_i], face), 1.0e-12, "outlet correction should add the configured flux trim"));
        PetscCall(PicurvAssertRealNear(1.0, ubcs[ubcs_k][ubcs_j][ubcs_i].x, 1.0e-12, "outlet handler should copy Ucat into Ubcs.x"));
        PetscCall(PicurvAssertRealNear(1.0, ubcs[ubcs_k][ubcs_j][ubcs_i].y, 1.0e-12, "outlet handler should copy Ucat into Ubcs.y"));
        PetscCall(PicurvAssertRealNear(1.0, ubcs[ubcs_k][ubcs_j][ubcs_i].z, 1.0e-12, "outlet handler should copy Ucat into Ubcs.z"));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(PicurvAssertRealNear(0.0, local_inflow, 1.0e-12, "outlet handler should not add to inflow"));
        PetscCall(PicurvAssertRealNear(1.2 * area, local_outflow, 1.0e-12, "outlet PostStep should report the corrected flux"));
 
        PetscCall(DestroyBoundaryHandler(&bc));
        ResetBoundaryFaceConfig(&user->boundary_faces[face]);
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests non-periodic cell-field finalization without modifying Ucont.
 */
static PetscErrorCode TestFinalizePostProjectionCellFieldsNonPeriodic(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Cmpnts ***ucat = NULL;
    Cmpnts ***lucat = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    PetscCall(VecSet(user->Ucat, 4.0));
    PetscCall(VecSet(user->Bcs.Ubcs, 1.0));
    PetscCall(VecSet(user->P, 3.0));
    PetscCall(VecSet(user->Ucont, 7.0));
 
    PetscCall(FinalizePostProjectionCellFields(user));
 
    PetscCall(PicurvAssertVecConstant(user->Ucont, 7.0, 1.0e-12,
                                      "post-projection cell finalization must preserve non-periodic Ucont"));
    PetscCall(PicurvAssertVecConstant(user->P, 3.0, 1.0e-12,
                                      "constant pressure should remain unchanged during cell finalization"));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lUcat, &lucat));
    PetscCall(PicurvAssertRealNear(-2.0, ucat[2][2][0].x, 1.0e-12,
                                    "non-periodic Ucat dummy face should use Ubcs extrapolation"));
    PetscCall(PicurvAssertRealNear(-2.0, lucat[2][2][0].x, 1.0e-12,
                                    "non-periodic local Ucat should contain the finalized dummy face"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lUcat, &lucat));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
/**
 * @brief EnforceRHSBoundaryConditions() must zero a row exactly when it is not an unknown.
 *
 * This is the invariant that ties the residual path to ClassifyMomentumRow(). It is
 * checked across a non-periodic, a singly-periodic and a fully-periodic configuration
 * because the defect it guards against was visible only when an axis was periodic:
 * the duplicate column at index 0 was left unzeroed, so ComputeTotalResidual()'s BDF
 * term accumulated there on every call and the residual norm stopped describing the
 * state. A sentinel fill makes an untouched row detectable rather than accidentally
 * correct.
 */
static PetscErrorCode TestEnforceRHSBoundaryConditionsMatchesRowClassification(void)
{
    const struct {
        const char *label;
        PetscBool   x_periodic, y_periodic, z_periodic;
    } configs[] = {
        {"non-periodic",     PETSC_FALSE, PETSC_FALSE, PETSC_FALSE},
        {"xi-periodic",      PETSC_TRUE,  PETSC_FALSE, PETSC_FALSE},
        {"fully-periodic",   PETSC_TRUE,  PETSC_TRUE,  PETSC_TRUE },
    };
    const PetscScalar sentinel = 7.0;
 
    PetscFunctionBeginUser;
    for (size_t c = 0; c < sizeof(configs) / sizeof(configs[0]); ++c) {
        SimCtx  *simCtx = NULL;
        UserCtx *user = NULL;
        Cmpnts ***rhs = NULL;
        const PetscBool periodic[3] = {configs[c].x_periodic, configs[c].y_periodic, configs[c].z_periodic};
        const BCFace neg_face[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
        const BCFace pos_face[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
 
        PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 6, 6, 6,
                                                             configs[c].x_periodic,
                                                             configs[c].y_periodic,
                                                             configs[c].z_periodic));
        for (PetscInt axis = 0; axis < 3; ++axis) {
            const BCType type = periodic[axis] ? PERIODIC : WALL;
            user->boundary_faces[neg_face[axis]].mathematical_type = type;
            user->boundary_faces[pos_face[axis]].mathematical_type = type;
        }
 
        PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
        PetscCall(VecSet(user->Rhs, sentinel));
        PetscCall(EnforceRHSBoundaryConditions(user));
 
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Rhs, &rhs));
        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 PetscScalar *row = &rhs[k][j][i].x;
                    for (PetscInt component = 0; component < 3; ++component) {
                        PetscInt ri, rj, rk;
                        const MomentumRowType type =
                            ClassifyMomentumRow(user, i, j, k, component, &ri, &rj, &rk);
                        const PetscScalar expected = (type == MOM_ROW_PHYSICAL) ? sentinel : 0.0;
                        char context[192];
 
                        PetscCall(PetscSNPrintf(context, sizeof(context),
                            "%s: row (i=%d,j=%d,k=%d,c=%d) classified %d must be %s",
                            configs[c].label, (int)i, (int)j, (int)k, (int)component, (int)type,
                            (type == MOM_ROW_PHYSICAL) ? "left untouched" : "zeroed"));
                        PetscCall(PicurvAssertRealNear((PetscReal)expected, (PetscReal)row[component],
                                                       1.0e-14, context));
                    }
                }
            }
        }
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Rhs, &rhs));
 
        PetscCall(VecDestroy(&user->Rhs));
        PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Allocates the wall-model storage the way production allocates it.
 *
 * Production creates each global on `da` and ghosts it into the local view. A fixture
 * that allocates differently - a different DM, only one of a pair, or only one of the
 * fields the same block creates - tests a layout the solver never uses, which is how a
 * DM/dof mismatch survived in `ApplyWallFunction` while these cases passed. The wall
 * eddy viscosity is allocated in that same block and cleared at the top of every pass,
 * so a fixture omitting it dereferences a null vector before reaching any assertion.
 * Every wall case routes through here so the two cannot drift.
 */
static PetscErrorCode SeedWallFrictionVelocityStorage(UserCtx *user)
{
    PetscFunctionBeginUser;
    PetscCall(DMCreateGlobalVector(user->da, &user->Friction_Velocity));
    PetscCall(VecSet(user->Friction_Velocity, 0.0));
    PetscCall(DMCreateLocalVector(user->da, &user->lFriction_Velocity));
    PetscCall(VecSet(user->lFriction_Velocity, 0.0));
    PetscCall(DMCreateGlobalVector(user->da, &user->Nu_Wall));
    PetscCall(VecSet(user->Nu_Wall, 0.0));
    PetscCall(DMCreateLocalVector(user->da, &user->lNu_Wall));
    PetscCall(VecSet(user->lNu_Wall, 0.0));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Seeds the fields ApplyWallFunction() reads on a single wall face.
 *
 * A uniform tangential velocity and unit metrics make the wall distances and the
 * face normal analytically known, so the test can reproduce the model's own
 * arithmetic without copying the integration's index choices.
 */
static PetscErrorCode SeedWallFunctionFixture(UserCtx *user, PetscReal tangential_speed)
{
    Cmpnts ***ucat = NULL, ***csi = NULL;
    PetscReal ***aj = NULL;
    DMDALocalInfo info = user->info;
 
    PetscFunctionBeginUser;
    PetscCall(VecSet(user->Aj, 1.0));
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecGetArray(user->fda, user->Csi, &csi));
    PetscCall(DMDAVecGetArray(user->da, user->Aj, &aj));
    for (PetscInt k = info.zs; k < info.zs + info.zm; ++k)
        for (PetscInt j = info.ys; j < info.ys + info.ym; ++j)
            for (PetscInt i = info.xs; i < info.xs + info.xm; ++i) {
                /* Flow along eta, so the whole velocity is tangential to a xi wall. */
                ucat[k][j][i] = (Cmpnts){.x = 0.0, .y = tangential_speed, .z = 0.0};
                csi[k][j][i] = (Cmpnts){.x = 1.0, .y = 0.0, .z = 0.0};
                aj[k][j][i] = 1.0;
            }
    PetscCall(DMDAVecRestoreArray(user->da, user->Aj, &aj));
    PetscCall(DMDAVecRestoreArray(user->fda, user->Csi, &csi));
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Aj, INSERT_VALUES, user->lAj));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Aj, INSERT_VALUES, user->lAj));
    PetscCall(DMGlobalToLocalBegin(user->fda, user->Csi, INSERT_VALUES, user->lCsi));
    PetscCall(DMGlobalToLocalEnd(user->fda, user->Csi, INSERT_VALUES, user->lCsi));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Tests that the wall-function integration reproduces the model it dispatches to.
 *
 * The kernels are covered separately; what is unverified without this is the wiring
 * around them, where the plausible failures are silent. Feeding the reference and
 * boundary wall distances in the wrong order, taking the reference velocity from the
 * wrong cell, or normalising the face normal incorrectly all still produce a
 * plausible-looking corrected velocity. The check here is that the corrected
 * first-interior velocity is exactly what the log-law model predicts at that cell's
 * distance, for the friction velocity the integration itself stored.
 */
static PetscErrorCode TestApplyWallFunctionMatchesTheLogLawItDispatches(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL, ***ubcs = NULL;
    PetscReal ***utau = NULL;
    const PetscReal speed = 1.0;
    const PetscReal roughness = 1.0e-16;   /* smooth wall */
    PetscReal    kinematic_viscosity, distance_first, distance_second;
    PetscReal    stored_utau, corrected_speed, expected_speed, reference_check;
    PetscReal    untouched_before = 0.0, untouched_after = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->ren = 1000.0;
    simCtx->wall_roughness_height = roughness;
    kinematic_viscosity = 1.0 / simCtx->ren;
 
    PetscCall(SeedWallFrictionVelocityStorage(user));
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = WALL;
    PetscCall(SeedWallFunctionFixture(user, speed));
 
    /* With unit Jacobians and a unit face area the model's distances are exact:
       the first cell centre sits half a cell from the wall, the second one and a
       half. Reproducing them here rather than reading them back is what makes a
       swapped-distance wiring error detectable. */
    distance_first = 0.5;
    distance_second = 2.0 * distance_first + 0.5;
 
    /* Disabled: the routine must not touch the field at all. */
    simCtx->wallfunction = 0;
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    untouched_before = ucat[2][2][1].y;
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    PetscCall(ApplyWallFunction(user));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    untouched_after = ucat[2][2][1].y;
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    PetscCall(PicurvAssertRealNear(untouched_before, untouched_after, 0.0,
                                   "a disabled wall function must leave the velocity untouched"));
 
    /* Enabled: the first interior cell is corrected onto the modelled profile. */
    simCtx->wallfunction = 1;
    PetscCall(ApplyWallFunction(user));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->Friction_Velocity, (const PetscReal ***)&utau));
    stored_utau = utau[2][2][1];
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Friction_Velocity, (const PetscReal ***)&utau));
    PetscCall(PicurvAssertBool((PetscBool)(stored_utau > 0.0),
                               "the integration must store a positive friction velocity"));
 
    /* The stored friction velocity must be the one that reproduces the *reference*
       cell's speed at the *reference* distance. Swapping the two distances passes the
       next check but fails this one. */
    reference_check = u_hydset_roughness(kinematic_viscosity, distance_second,
                                         stored_utau, roughness);
    PetscCall(PicurvAssertRealNear(speed, reference_check, 1.0e-6,
                                   "the friction velocity must invert the log law at the "
                                   "reference cell's distance"));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    corrected_speed = ucat[2][2][1].y;
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
 
    expected_speed = u_hydset_roughness(kinematic_viscosity, distance_first,
                                        stored_utau, roughness);
    PetscCall(PicurvAssertRealNear(expected_speed, corrected_speed, 1.0e-9,
                                   "the corrected velocity must be the modelled profile at "
                                   "the first interior cell's distance"));
    PetscCall(PicurvAssertBool((PetscBool)(corrected_speed < speed),
                               "a log-law wall correction must slow the first interior cell "
                               "relative to the reference cell"));
 
    /* The ghost row carries no boundary velocity of its own. */
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, (const Cmpnts ***)&ubcs));
    PetscCall(PicurvAssertRealNear(0.0, ubcs[2][2][0].y, 0.0,
                                   "the wall ghost row must carry a zero boundary velocity"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, (const Cmpnts ***)&ubcs));
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the model selector actually changes which law is applied.
 *
 * The selector was previously validated in the configuration layer and then dropped, so
 * every case ran the log law whatever it asked for. Distinct laws must produce distinct
 * corrections from identical inputs, or the selector is decorative again.
 */
static PetscErrorCode TestWallFunctionModelSelectorChoosesTheLaw(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal    corrected[4] = {0.0, 0.0, 0.0, 0.0};
    const PetscInt models[3] = {WALL_FUNCTION_LOG_LAW, WALL_FUNCTION_WERNER, WALL_FUNCTION_CABOT};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->ren = 1000.0;
    simCtx->wall_roughness_height = 1.0e-16;
    PetscCall(SeedWallFrictionVelocityStorage(user));
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = WALL;
 
    for (PetscInt m = 0; m < 3; ++m) {
        PetscCall(SeedWallFunctionFixture(user, 1.0));
        simCtx->wallfunction = models[m];
        PetscCall(ApplyWallFunction(user));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
        corrected[m] = ucat[2][2][1].y;
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
        PetscCall(PicurvAssertBool((PetscBool)(corrected[m] > 0.0 && corrected[m] < 1.0),
                                   "every wall model must slow the first interior cell"));
    }
 
    /* Each law has its own profile, so no two may coincide. */
    PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(corrected[0] - corrected[1]) > 1.0e-9),
                               "the Werner-Wengle model must differ from the log law"));
    PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(corrected[0] - corrected[2]) > 1.0e-9),
                               "the Cabot model must differ from the log law"));
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that a face which is not a wall is left alone by the wall function.
 */
static PetscErrorCode TestApplyWallFunctionSkipsNonWallFaces(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    Cmpnts    ***ucat = NULL;
    PetscReal    before = 0.0, after = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->ren = 1000.0;
    simCtx->wallfunction = 1;
    simCtx->wall_roughness_height = 1.0e-16;
    PetscCall(SeedWallFrictionVelocityStorage(user));
    /* Every face left at its default, so none is a wall. */
    PetscCall(SeedWallFunctionFixture(user, 1.0));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    before = ucat[2][2][1].y;
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    PetscCall(ApplyWallFunction(user));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
    after = ucat[2][2][1].y;
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, (const Cmpnts ***)&ucat));
 
    PetscCall(PicurvAssertRealNear(before, after, 0.0,
                                   "a face that is not a wall must not be corrected"));
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that every wall law's friction velocity responds to the flow it is given.
 *
 * A law that ignores the resolved velocity and reports a fixed number still slows the
 * first cell and still differs from the other laws, so the selector test above passes
 * over it. Driving a law at two different speeds and requiring the friction velocity to
 * follow is what distinguishes a model from a constant.
 *
 * All three laws are covered. Werner-Wengle is the reason this case exists: it used to
 * report its Newton initial guess rather than solving, so every wall cell in a run
 * carried the same u_tau while still slowing the cell and still differing from the other
 * laws, which is everything the selector test above checks.
 */
static PetscErrorCode TestWallModelFrictionVelocityRespondsToTheFlow(void)
{
    SimCtx        *simCtx = NULL;
    UserCtx       *user = NULL;
    PetscReal   ***utau = NULL;
    const PetscInt models[3] = {WALL_FUNCTION_LOG_LAW, WALL_FUNCTION_WERNER,
                                WALL_FUNCTION_CABOT};
    const char    *names[3] = {"the log law must respond to the flow it is given",
                               "Werner-Wengle must respond to the flow it is given",
                               "Cabot must respond to the flow it is given"};
    const PetscReal speeds[2] = {0.5, 2.0};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->ren = 1000.0;
    simCtx->wall_roughness_height = 1.0e-16;
    PetscCall(SeedWallFrictionVelocityStorage(user));
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = WALL;
 
    for (PetscInt m = 0; m < 3; ++m) {
        PetscReal stored[2] = {0.0, 0.0};
 
        simCtx->wallfunction = models[m];
        for (PetscInt s = 0; s < 2; ++s) {
            PetscCall(SeedWallFunctionFixture(user, speeds[s]));
            PetscCall(ApplyWallFunction(user));
            PetscCall(DMDAVecGetArrayRead(user->da, user->Friction_Velocity,
                                          (const PetscReal ***)&utau));
            stored[s] = utau[2][2][1];
            PetscCall(DMDAVecRestoreArrayRead(user->da, user->Friction_Velocity,
                                              (const PetscReal ***)&utau));
            PetscCall(PicurvAssertBool((PetscBool)(stored[s] > 0.0),
                                       "a wall law must report a positive friction velocity"));
        }
 
        /* Faster flow, more wall stress. Any law that reports its initial guess, or any
           other constant, fails here whatever that constant is. */
        PetscCall(PicurvAssertBool((PetscBool)(stored[1] > stored[0] * 1.05),
                                   names[m]));
    }
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the friction-velocity storage matches what the field catalog declares.
 *
 * The wall model reaches its storage through a DM, and PETSc only rejects a mismatched
 * DM/vector pairing at the moment of access - inside a branch that runs only when wall
 * functions are enabled on a case that has a WALL face. That is how storage created from
 * a dof-3 DM and read through a dof-1 one stayed latent. Comparing the allocation against
 * the catalog descriptor turns that into a failure here.
 */
static PetscErrorCode TestWallFrictionVelocityStorageMatchesTheCatalog(void)
{
    SimCtx                *simCtx = NULL;
    UserCtx               *user = NULL;
    const FieldDescriptor *descriptor = NULL;
    FieldView              view;
    PetscInt               global_size = 0, local_size = 0;
    PetscInt               xm = 0, ym = 0, zm = 0, gxm = 0, gym = 0, gzm = 0;
    PetscReal           ***utau = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    PetscCall(SeedWallFrictionVelocityStorage(user));
 
    PetscCall(FieldGetDescriptor(FIELD_ID_U_TAU, &descriptor));
    PetscCall(PicurvAssertBool((PetscBool)(descriptor->dof == 1),
                               "the friction velocity is one scalar per cell"));
    PetscCall(PicurvAssertBool((PetscBool)(descriptor->dm_kind == FIELD_DM_DA),
                               "a scalar cell field belongs on the dof-1 DM"));
 
    PetscCall(FieldGetView(user, FIELD_ID_U_TAU, &view));
    PetscCall(PicurvAssertBool((PetscBool)(view.global_vec == user->Friction_Velocity &&
                                           view.local_vec == user->lFriction_Velocity),
                               "the catalog must resolve to the vectors setup allocated"));
 
    /* The sizes the DM implies, against the sizes the vectors actually have. */
    PetscCall(DMDAGetCorners(user->da, NULL, NULL, NULL, &xm, &ym, &zm));
    PetscCall(DMDAGetGhostCorners(user->da, NULL, NULL, NULL, &gxm, &gym, &gzm));
    PetscCall(VecGetLocalSize(user->Friction_Velocity, &global_size));
    PetscCall(VecGetLocalSize(user->lFriction_Velocity, &local_size));
    PetscCall(PicurvAssertBool((PetscBool)(global_size == xm * ym * zm * descriptor->dof),
                               "the global friction velocity must match the owned extent"));
    PetscCall(PicurvAssertBool((PetscBool)(local_size == gxm * gym * gzm * descriptor->dof),
                               "the local friction velocity must match the ghosted extent"));
 
    /* The access the wall model itself performs. A mismatched pairing fails right here. */
    PetscCall(DMDAVecGetArray(user->da, user->Friction_Velocity, &utau));
    PetscCall(DMDAVecRestoreArray(user->da, user->Friction_Velocity, &utau));
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the wall pass records the statistics its diagnostics report.
 *
 * `y+` cannot be recovered downstream, so it is formed where the wall distance is still
 * in hand. This checks that what is accumulated there is the same quantity the reported
 * column claims to be, and that a pass over no wall face records nothing rather than
 * recording zeros.
 */
static PetscErrorCode TestWallModelDiagnosticsRecordTheCorrectedCells(void)
{
    SimCtx      *simCtx = NULL;
    UserCtx     *user = NULL;
    PetscReal    kinematic_viscosity = 0.0;
    PetscReal    expected_y_plus = 0.0;
    const PetscReal distance_first = 0.5;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
    simCtx->ren = 1000.0;
    simCtx->wall_roughness_height = 1.0e-16;
    kinematic_viscosity = 1.0 / simCtx->ren;
    PetscCall(SeedWallFrictionVelocityStorage(user));
    user->boundary_faces[BC_FACE_NEG_X].mathematical_type = WALL;
    PetscCall(SeedWallFunctionFixture(user, 1.0));
 
    /* Disabled: nothing is corrected, so nothing is recorded. */
    simCtx->wallfunction = 0;
    PetscCall(ApplyWallFunction(user));
    PetscCall(PicurvAssertBool((PetscBool)(user->wall_diagnostics.cells == 0),
                               "a disabled wall model must record no corrected cells"));
 
    simCtx->wallfunction = 1;
    PetscCall(ApplyWallFunction(user));
    PetscCall(PicurvAssertBool((PetscBool)(user->wall_diagnostics.cells > 0),
                               "an enabled wall model must record the cells it corrected"));
 
    /* One face, so the sums are one cell's values and y+ is exactly u_tau * y / nu.
       Reproducing it from the stored friction velocity is what makes a wrong distance
       or a wrong viscosity in the accumulation visible. */
    expected_y_plus = user->wall_diagnostics.friction_velocity_sum * distance_first /
                      kinematic_viscosity / (PetscReal)user->wall_diagnostics.cells;
    PetscCall(PicurvAssertRealNear(user->wall_diagnostics.y_plus_sum /
                                       (PetscReal)user->wall_diagnostics.cells,
                                   expected_y_plus, 1.0e-9,
                                   "y+ must be formed from the same distance the "
                                   "correction used"));
    PetscCall(PicurvAssertRealNear(user->wall_diagnostics.wall_distance_sum /
                                       (PetscReal)user->wall_diagnostics.cells,
                                   distance_first, 1.0e-9,
                                   "the recorded wall distance must be the first cell's"));
    PetscCall(PicurvAssertBool((PetscBool)(user->wall_diagnostics.friction_velocity_min <=
                                           user->wall_diagnostics.friction_velocity_max),
                               "the recorded extrema must bracket each other"));
 
    /* A second pass must describe that pass alone, not both. */
    {
        const PetscInt first_pass_cells = user->wall_diagnostics.cells;
 
        PetscCall(ApplyWallFunction(user));
        PetscCall(PicurvAssertBool((PetscBool)(user->wall_diagnostics.cells == first_pass_cells),
                                   "each pass must replace the previous pass's sample set"));
    }
 
    PetscCall(VecDestroy(&user->Friction_Velocity));
    PetscCall(VecDestroy(&user->lFriction_Velocity));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Checks each cell-centred field's dummy-cell rule, faces then edges.
 * @details Interior cells carry index-unique sentinels and every dummy cell starts
 *          poisoned, so a wrong offset or a missed face shows. Velocity dummies must put
 *          the face average on `Ubcs`; pressure and the particle scalar copy the adjacent
 *          cell (zero normal gradient); each edge dummy is the mean of its two face
 *          neighbours; and a field with no rule, or one that is not cell-centred, is refused.
 */
static PetscErrorCode TestDummyCellsFollowEachFieldsBoundaryRule(void)
{
    SimCtx       *simCtx = NULL;
    UserCtx      *user = NULL;
    const FieldId scalars[] = {FIELD_ID_P, FIELD_ID_PSI};
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 6, 7));
    const DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscCall(VecSet(user->Bcs.Ubcs, 0.0));
    {
        PetscReal ***p = NULL, ***psi = NULL;
        Cmpnts    ***ucat = NULL, ***ubcs = NULL;
 
        PetscCall(DMDAVecGetArray(user->da, user->P, &p));
        PetscCall(DMDAVecGetArray(user->da, user->Psi, &psi));
        PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
        PetscCall(DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs));
        for (PetscInt k = info.zs; k < info.zs + info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++)
        for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
            const PetscBool interior = (PetscBool)(i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1);
            const PetscReal sentinel = 100.0 * k + 10.0 * j + i;
 
            p[k][j][i] = interior ? sentinel : -999.0;
            psi[k][j][i] = interior ? -sentinel : -999.0;
            ucat[k][j][i].x = interior ? sentinel : -999.0;
            ucat[k][j][i].y = interior ? 2.0 * sentinel : -999.0;
            ucat[k][j][i].z = interior ? 3.0 * sentinel : -999.0;
            ubcs[k][j][i] = (Cmpnts){1.0, 2.0, 3.0};
        }
        PetscCall(DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs));
        PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
        PetscCall(DMDAVecRestoreArray(user->da, user->Psi, &psi));
        PetscCall(DMDAVecRestoreArray(user->da, user->P, &p));
    }
 
    PetscCall(UpdateDummyCells(user, FIELD_ID_UCAT));
    PetscCall(UpdateCornerNodes(user, FIELD_ID_UCAT));
    for (size_t n = 0; n < sizeof(scalars) / sizeof(scalars[0]); ++n) {
        PetscCall(UpdateDummyCells(user, scalars[n]));
        PetscCall(UpdateCornerNodes(user, scalars[n]));
    }
 
    for (size_t n = 0; n < sizeof(scalars) / sizeof(scalars[0]); ++n) {
        FieldView  view;
        PetscReal ***f = NULL;
 
        PetscCall(FieldGetView(user, scalars[n], &view));
        PetscCall(DMDAVecGetArrayRead(user->da, view.global_vec, &f));
        for (PetscInt k = 1; k < mz - 1; k++) for (PetscInt j = 1; j < my - 1; j++) {
            PetscCall(PicurvAssertRealNear(f[k][j][1], f[k][j][0], 0.0, "-x dummy copies the adjacent cell"));
            PetscCall(PicurvAssertRealNear(f[k][j][mx - 2], f[k][j][mx - 1], 0.0, "+x dummy copies the adjacent cell"));
        }
        for (PetscInt k = 1; k < mz - 1; k++) for (PetscInt i = 1; i < mx - 1; i++) {
            PetscCall(PicurvAssertRealNear(f[k][1][i], f[k][0][i], 0.0, "-y dummy copies the adjacent cell"));
            PetscCall(PicurvAssertRealNear(f[k][my - 2][i], f[k][my - 1][i], 0.0, "+y dummy copies the adjacent cell"));
        }
        for (PetscInt j = 1; j < my - 1; j++) for (PetscInt i = 1; i < mx - 1; i++) {
            PetscCall(PicurvAssertRealNear(f[1][j][i], f[0][j][i], 0.0, "-z dummy copies the adjacent cell"));
            PetscCall(PicurvAssertRealNear(f[mz - 2][j][i], f[mz - 1][j][i], 0.0, "+z dummy copies the adjacent cell"));
        }
        for (PetscInt j = 0; j < my; j++) {
            PetscCall(PicurvAssertRealNear(0.5 * (f[1][j][0] + f[0][j][1]), f[0][j][0], 0.0,
                                           "a (-z,-x) edge dummy is the mean of its face neighbours"));
        }
        PetscCall(PicurvAssertBool((PetscBool)(f[0][0][0] != -999.0 && f[mz - 1][my - 1][mx - 1] != -999.0),
                                   "corner dummies are written"));
        PetscCall(DMDAVecRestoreArrayRead(user->da, view.global_vec, &f));
    }
    {
        const Cmpnts ***ucat = NULL;
 
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
        for (PetscInt k = 1; k < mz - 1; k++) for (PetscInt j = 1; j < my - 1; j++) {
            PetscCall(PicurvAssertRealNear(1.0, 0.5 * (ucat[k][j][0].x + ucat[k][j][1].x), 1.0e-12,
                                           "a velocity dummy puts the face average on Ubcs"));
            PetscCall(PicurvAssertRealNear(3.0, 0.5 * (ucat[k][j][mx - 1].z + ucat[k][j][mx - 2].z), 1.0e-12,
                                           "a velocity dummy puts the face average on Ubcs"));
        }
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
    }
 
    {
        PetscErrorCode refused;
 
        PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
        refused = UpdateDummyCells(user, FIELD_ID_NVERT);
        PetscCall(PetscPopErrorHandler());
        PetscCall(PicurvAssertIntEqual(PETSC_ERR_SUP, refused, "a field with no boundary rule is refused"));
        PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
        refused = UpdateDummyCells(user, FIELD_ID_COORDINATES);
        PetscCall(PetscPopErrorHandler());
        PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONG, refused, "a field that is not cell-centred is refused"));
    }
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Runs the unit-boundaries PETSc test binary.
 */
 
int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"can-rank-service-face-matches-inlet-when-defined", TestCanRankServiceFaceMatchesInletWhenDefined},
        {"can-rank-service-inlet-face-requires-definition", TestCanRankServiceInletFaceRequiresDefinition},
        {"boundary-condition-factory-assignments", TestBoundaryConditionFactoryAssignments},
        {"boundary-condition-factory-implemented-handler-matrix", TestBoundaryConditionFactoryImplementedHandlerMatrix},
        {"boundary-condition-factory-rejects-unsupported-handler", TestBoundaryConditionFactoryRejectsUnsupportedHandler},
        {"deterministic-face-grid-location-matrix", TestGetDeterministicFaceGridLocationFaceMatrix},
        {"random-inlet-face-location-matrix", TestGetRandomCellAndLogicalCoordsOnInletFaceMatrix},
        {"wall-no-slip-handler-face-matrix", TestWallNoSlipHandlerFaceMatrix},
        {"dummy-cells-follow-each-fields-boundary-rule", TestDummyCellsFollowEachFieldsBoundaryRule},
        {"inlet-constant-velocity-handler-behavior", TestInletConstantVelocityHandlerBehavior},
        {"inlet-common-application-curvilinear-metric", TestInletCommonApplicationCurvilinearMetric},
        {"inlet-constant-velocity-handler-face-matrix", TestInletConstantVelocityHandlerFaceMatrix},
        {"inlet-parabolic-profile-handler-behavior", TestInletParabolicProfileHandlerBehavior},
        {"inlet-parabolic-profile-handler-face-matrix", TestInletParabolicProfileHandlerFaceMatrix},
        {"inlet-parabolic-profile-uniform-along-periodic-axis", TestInletParabolicProfileUniformAlongPeriodicAxis},
        {"inlet-profile-from-file-handler-behavior", TestInletProfileFromFileHandlerBehavior},
        {"outlet-conservation-handler-behavior", TestOutletConservationHandlerBehavior},
        {"outlet-conservation-handler-face-matrix", TestOutletConservationHandlerFaceMatrix},
        {"finalize-post-projection-cell-fields-non-periodic", TestFinalizePostProjectionCellFieldsNonPeriodic},
        {"enforce-rhs-bcs-matches-row-classification", TestEnforceRHSBoundaryConditionsMatchesRowClassification},
        {"apply-wall-function-matches-the-log-law", TestApplyWallFunctionMatchesTheLogLawItDispatches},
        {"apply-wall-function-skips-non-wall-faces", TestApplyWallFunctionSkipsNonWallFaces},
        {"wall-function-model-selector-chooses-the-law", TestWallFunctionModelSelectorChoosesTheLaw},
        {"wall-model-friction-velocity-responds-to-the-flow", TestWallModelFrictionVelocityRespondsToTheFlow},
        {"wall-friction-velocity-storage-matches-the-catalog", TestWallFrictionVelocityStorageMatchesTheCatalog},
        {"wall-model-diagnostics-record-the-corrected-cells", TestWallModelDiagnosticsRecordTheCorrectedCells},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv boundary tests");
    if (ierr) {
        return (int)ierr;
    }
 
    ierr = PicurvRunTests("unit-boundaries", cases, sizeof(cases) / sizeof(cases[0]));
    if (ierr) {
        PetscFinalize();
        return (int)ierr;
    }
 
    ierr = PetscFinalize();
    return (int)ierr;
}