test_momentum_newton_krylov.c

Go to the documentation of this file.

/**
 * @file test_momentum_newton_krylov.c
 * @brief Focused tests for the version-one matrix-free momentum solver.
 *
 * The implementation is included intentionally: its callback helpers remain
 * private in production while this translation unit can verify them directly.
 */
 
#include "test_support.h"
#include "initialcondition.h"
 
/* Rename only the included public entry points. Private callback tests use this
 * copy, while lifecycle tests below call the separately linked production object. */
#define MomentumSolver_NewtonKrylov MomentumSolver_NewtonKrylov_PrivateCopy
#define MomentumSolver_NewtonKrylov_Destroy MomentumSolver_NewtonKrylov_Destroy_PrivateCopy
#include "../../src/momentum_newton_krylov.c"
#undef MomentumSolver_NewtonKrylov
#undef MomentumSolver_NewtonKrylov_Destroy
 
static const char *geometric_periodic_bcs =
    "-Xi PERIODIC geometric\n"
    "+Xi PERIODIC geometric\n"
    "-Eta WALL noslip\n"
    "+Eta WALL noslip\n"
    "-Zeta INLET prescribed_velocity source_id=0\n"
    "+Zeta OUTLET conservation\n";
 
/* Two free-slip walls facing two no-slip ones: the free-slip boundary value is read
 * from the trial velocity, so the residual must still be a function of X alone. */
static const char *free_slip_wall_bcs =
    "-Xi WALL free_slip\n"
    "+Xi WALL free_slip\n"
    "-Eta WALL noslip\n"
    "+Eta WALL noslip\n"
    "-Zeta WALL noslip\n"
    "+Zeta WALL noslip\n";
 
static const char *fixed_wall_bcs =
    "-Xi WALL noslip\n"
    "+Xi WALL noslip\n"
    "-Eta WALL noslip\n"
    "+Eta WALL noslip\n"
    "-Zeta WALL noslip\n"
    "+Zeta WALL noslip\n";
 
static const char *parabolic_bcs =
    "-Xi WALL noslip\n"
    "+Xi WALL noslip\n"
    "-Eta WALL noslip\n"
    "+Eta WALL noslip\n"
    "-Zeta INLET prescribed_velocity source_id=1\n"
    "+Zeta OUTLET conservation\n";
 
static const char *periodic_x_bcs =
    "-Xi PERIODIC geometric\n+Xi PERIODIC geometric\n"
    "-Eta WALL noslip\n+Eta WALL noslip\n-Zeta WALL noslip\n+Zeta WALL noslip\n";
 
static const char *periodic_y_bcs =
    "-Xi WALL noslip\n+Xi WALL noslip\n"
    "-Eta PERIODIC geometric\n+Eta PERIODIC geometric\n-Zeta WALL noslip\n+Zeta WALL noslip\n";
 
static const char *periodic_z_bcs =
    "-Xi WALL noslip\n+Xi WALL noslip\n-Eta WALL noslip\n+Eta WALL noslip\n"
    "-Zeta PERIODIC geometric\n+Zeta PERIODIC geometric\n";
 
static const char *periodic_xy_bcs =
    "-Xi PERIODIC geometric\n+Xi PERIODIC geometric\n"
    "-Eta PERIODIC geometric\n+Eta PERIODIC geometric\n"
    "-Zeta WALL noslip\n+Zeta WALL noslip\n";
 
static const char *periodic_xyz_bcs =
    "-Xi PERIODIC geometric\n+Xi PERIODIC geometric\n"
    "-Eta PERIODIC geometric\n+Eta PERIODIC geometric\n"
    "-Zeta PERIODIC geometric\n+Zeta PERIODIC geometric\n";
 
/** @brief Checks a structured log's row count and required text after a collective solve. */
static PetscErrorCode AssertNewtonLog(const char *path, PetscInt expected_rows,
                                      const char *needle_a, const char *needle_b)
{
    FILE *file = NULL;
    char line[4096];
    PetscInt rows = 0;
    PetscBool found_a = needle_a ? PETSC_FALSE : PETSC_TRUE;
    PetscBool found_b = needle_b ? PETSC_FALSE : PETSC_TRUE;
 
    PetscFunctionBeginUser;
    PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
    file = fopen(path, "r");
    PetscCheck(file != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
               "Expected Newton log does not exist: %s", path);
    while (fgets(line, sizeof(line), file)) {
        if (line[0] >= '0' && line[0] <= '9') rows++; /* CSV data rows; the header and comments are not */
        if (needle_a && strstr(line, needle_a)) found_a = PETSC_TRUE;
        if (needle_b && strstr(line, needle_b)) found_b = PETSC_TRUE;
    }
    fclose(file);
    if (expected_rows >= 0) {
        PetscCall(PicurvAssertIntEqual(expected_rows, rows,
                                       "Newton log must contain one nonduplicated row per solve"));
    } else {
        PetscCall(PicurvAssertBool((PetscBool)(rows >= -expected_rows),
                                   "enabled Newton history must contain the expected iteration rows"));
    }
    if (!found_a || !found_b) {
        PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s lacks '%s'\n", path, !found_a ? needle_a : needle_b));
    }
    PetscCall(PicurvAssertBool(found_a, "Newton log is missing required structured content"));
    PetscCall(PicurvAssertBool(found_b, "Newton log is missing required structured content"));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Builds and initializes a small runtime context for Newton tests.
 * @param bcs Optional boundary configuration text.
 * @param simCtx Returned simulation context.
 * @param user Returned finest-level block context.
 * @param tmpdir Returned temporary directory.
 * @param tmpdir_len Capacity of tmpdir.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode BuildNewtonFixture(const char *bcs, SimCtx **simCtx, UserCtx **user,
                                         char *tmpdir, size_t tmpdir_len)
{
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContext(bcs, PETSC_FALSE, simCtx, user, tmpdir, tmpdir_len));
    PetscCall(InitializeEulerianState(*simCtx));
    (*simCtx)->mom_solver_type = MOMENTUM_SOLVER_NEWTON_KRYLOV;
    PetscCall(PicurvAssertBool((PetscBool)((*user)->Rhs == NULL),
                               "Newton fixture must enter with no persistent Rhs workspace"));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Destroys a Newton test fixture and its temporary files.
 * @param simCtx Fixture simulation context.
 * @param tmpdir Fixture temporary directory.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode DestroyNewtonFixture(SimCtx **simCtx, char *tmpdir)
{
    PetscFunctionBeginUser;
    PetscCall(PicurvDestroyRuntimeContext(simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Writes one static 5x5 PICSLICE profile used by the full runtime fixture.
 * @param path Output profile path.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode WriteNewtonPicSlice(const char *path)
{
    FILE *fd = NULL;
 
    PetscFunctionBeginUser;
    fd = fopen(path, "w");
    PetscCheck(fd != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
               "Could not create Newton test PICSLICE %s.", path);
    PetscCheck(fprintf(fd, "PICSLICE\n1\n5 5\n") >= 0,
               PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE, "Could not write PICSLICE header.");
    for (PetscInt row = 0; row < 25; ++row) {
        PetscCheck(fprintf(fd, "%.16e\n", 1.0 + 0.01 * (double)row) >= 0,
                   PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE, "Could not write PICSLICE value.");
    }
    fclose(fd);
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Checks callback repeatability, diagnostic-state independence, and X integrity.
 * @param bcs Boundary configuration text, or NULL for the standard inlet/outlet fixture.
 * @param label Configuration label used in assertion diagnostics.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode CheckResidualRepeatabilityForBC(const char *bcs, const char *label)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL, x_copy = NULL, f1 = NULL, f2 = NULL, delta = NULL;
    PetscReal norm = 0.0;
    MomentumNewtonKrylovContext ctx;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &x_copy));
    PetscCall(VecDuplicate(user->Ucont, &f1));
    PetscCall(VecDuplicate(user->Ucont, &f2));
    PetscCall(VecDuplicate(user->Ucont, &delta));
    PetscCall(VecCopy(user->Ucont, x));
    PetscCall(VecShift(x, 0.125));
    PetscCall(VecCopy(x, x_copy));
    ctx.user = user;
 
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f1, &ctx));
    /* Poison every piece of hidden state a non-deterministic residual could lean
     * on: the boundary flux/area diagnostics AND the persistent Cartesian fields
     * (Ucat/lUcat). The conservation-outlet handler reads lUcat during its first
     * boundary sweep, so a callback that does not reconstruct the Cartesian state
     * from X before applying boundary conditions would produce a different F here.
     * This assertion therefore fails if the deterministic pre-boundary seed in
     * MomentumNewtonKrylov_FormResidual() is ever removed. */
    simCtx->FluxInSum = 1234.0;
    simCtx->FluxOutSum = -4321.0;
    simCtx->FarFluxInSum = 77.0;
    simCtx->FarFluxOutSum = -88.0;
    PetscCall(VecSet(user->Ucat, 7.0));
    PetscCall(VecSet(user->lUcat, 7.0));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f2, &ctx));
    PetscCall(VecWAXPY(delta, -1.0, f1, f2));
    PetscCall(VecNorm(delta, NORM_INFINITY, &norm));
    PetscCheck(norm <= 1.0e-13, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "%s residual changed between identical evaluations (inf norm=%g).", label, (double)norm);
    PetscCall(VecNorm(delta, NORM_2, &norm));
    PetscCheck(norm <= 1.0e-12, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "%s residual changed between identical evaluations (L2 norm=%g).", label, (double)norm);
    PetscCall(VecWAXPY(delta, -1.0, x_copy, x));
    PetscCall(VecNorm(delta, NORM_INFINITY, &norm));
    PetscCheck(norm == 0.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "%s residual callback modified X (norm=%g).", label, (double)norm);
 
    PetscCall(VecDestroy(&delta));
    PetscCall(VecDestroy(&f2));
    PetscCall(VecDestroy(&f1));
    PetscCall(VecDestroy(&x_copy));
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Returns whether this rank owns one global DMDA grid point. */
static PetscBool OwnsStoredPoint(UserCtx *user, PetscInt i, PetscInt j, PetscInt k)
{
    return (PetscBool)(i >= user->info.xs && i < user->info.xs + user->info.xm &&
                       j >= user->info.ys && j < user->info.ys + user->info.ym &&
                       k >= user->info.zs && k < user->info.zs + user->info.zm);
}
 
/**
 * @brief Adds a scalar perturbation to one stored staggered component.
 * @param user Block context defining vector ownership.
 * @param vec Vector to modify.
 * @param i Global i index.
 * @param j Global j index.
 * @param k Global k index.
 * @param component Component 0=x, 1=y, 2=z.
 * @param delta Increment to apply.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode PerturbStoredValue(UserCtx *user, Vec vec, PetscInt i, PetscInt j,
                                         PetscInt k, PetscInt component, PetscScalar delta)
{
    Cmpnts ***a = NULL;
 
    PetscFunctionBeginUser;
    if (OwnsStoredPoint(user, i, j, k)) {
        PetscCall(DMDAVecGetArray(user->fda, vec, &a));
        if (component == 0) a[k][j][i].x += delta;
        else if (component == 1) a[k][j][i].y += delta;
        else a[k][j][i].z += delta;
        PetscCall(DMDAVecRestoreArray(user->fda, vec, &a));
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Reads one globally indexed stored component on any MPI decomposition.
 * @param user Block context defining vector ownership.
 * @param vec Vector to inspect.
 * @param i Global i index.
 * @param j Global j index.
 * @param k Global k index.
 * @param component Component 0=x, 1=y, 2=z.
 * @param value Returned globally reduced scalar.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode GetStoredValue(UserCtx *user, Vec vec, PetscInt i, PetscInt j,
                                     PetscInt k, PetscInt component, PetscScalar *value)
{
    Cmpnts ***a = NULL;
    PetscScalar local = 0.0;
 
    PetscFunctionBeginUser;
    if (OwnsStoredPoint(user, i, j, k)) {
        PetscCall(DMDAVecGetArrayRead(user->fda, vec, &a));
        local = component == 0 ? a[k][j][i].x : (component == 1 ? a[k][j][i].y : a[k][j][i].z);
        PetscCall(DMDAVecRestoreArrayRead(user->fda, vec, &a));
    }
    PetscCallMPI(MPI_Allreduce(&local, value, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Finite-differences one callback row with respect to one stored unknown.
 * @param user Active block context with allocated Rhs.
 * @param x Base trial vector.
 * @param row_i Residual-row i index.
 * @param row_j Residual-row j index.
 * @param row_k Residual-row k index.
 * @param row_component Residual-row component.
 * @param col_i Perturbed unknown i index.
 * @param col_j Perturbed unknown j index.
 * @param col_k Perturbed unknown k index.
 * @param col_component Perturbed unknown component.
 * @param derivative Returned finite-difference derivative of the row w.r.t. the unknown.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode MeasureStoredDerivative(UserCtx *user, Vec x,
                                              PetscInt row_i, PetscInt row_j, PetscInt row_k,
                                              PetscInt row_component, PetscInt col_i, PetscInt col_j,
                                              PetscInt col_k, PetscInt col_component,
                                              PetscReal *derivative)
{
    const PetscReal epsilon = 1.0e-6;
    Vec f0 = NULL, fp = NULL, xp = NULL;
    PetscScalar base_value = 0.0, perturbed_value = 0.0;
    MomentumNewtonKrylovContext ctx = {user};
 
    PetscFunctionBeginUser;
    PetscCall(VecDuplicate(x, &f0));
    PetscCall(VecDuplicate(x, &fp));
    PetscCall(VecDuplicate(x, &xp));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    PetscCall(VecCopy(x, xp));
    PetscCall(PerturbStoredValue(user, xp, col_i, col_j, col_k, col_component, epsilon));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, xp, fp, &ctx));
    PetscCall(GetStoredValue(user, f0, row_i, row_j, row_k, row_component, &base_value));
    PetscCall(GetStoredValue(user, fp, row_i, row_j, row_k, row_component, &perturbed_value));
    *derivative = PetscRealPart((perturbed_value - base_value) / epsilon);
    PetscCall(VecDestroy(&xp));
    PetscCall(VecDestroy(&fp));
    PetscCall(VecDestroy(&f0));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Asserts one finite-differenced callback row derivative equals an expected value.
 * @param user Active block context with allocated Rhs.
 * @param x Base trial vector.
 * @param row_i Residual-row i index.
 * @param row_j Residual-row j index.
 * @param row_k Residual-row k index.
 * @param row_component Residual-row component.
 * @param col_i Perturbed unknown i index.
 * @param col_j Perturbed unknown j index.
 * @param col_k Perturbed unknown k index.
 * @param col_component Perturbed unknown component.
 * @param expected Expected derivative.
 * @param tolerance Absolute derivative tolerance.
 * @param label Assertion label.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode CheckStoredDerivative(UserCtx *user, Vec x,
                                            PetscInt row_i, PetscInt row_j, PetscInt row_k,
                                            PetscInt row_component, PetscInt col_i, PetscInt col_j,
                                            PetscInt col_k, PetscInt col_component,
                                            PetscReal expected, PetscReal tolerance, const char *label)
{
    PetscReal derivative = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(MeasureStoredDerivative(user, x, row_i, row_j, row_k, row_component,
                                      col_i, col_j, col_k, col_component, &derivative));
    PetscCall(PicurvAssertRealNear(expected, derivative, tolerance, label));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies repeatable callback output and read-only trial input. */
static PetscErrorCode TestResidualRepeatabilityAndInputIntegrity(void)
{
    char profile_dir[PETSC_MAX_PATH_LEN] = "";
    char profile_path[PETSC_MAX_PATH_LEN];
    char file_bcs[2 * PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PetscCall(CheckResidualRepeatabilityForBC(fixed_wall_bcs, "fixed walls"));
    PetscCall(CheckResidualRepeatabilityForBC(free_slip_wall_bcs, "free-slip walls"));
    PetscCall(CheckResidualRepeatabilityForBC(NULL, "constant inlet/conservation outlet"));
    PetscCall(CheckResidualRepeatabilityForBC(parabolic_bcs, "parabolic inlet/conservation outlet"));
    PetscCall(CheckResidualRepeatabilityForBC(periodic_x_bcs, "x periodic"));
    PetscCall(CheckResidualRepeatabilityForBC(periodic_y_bcs, "y periodic"));
    PetscCall(CheckResidualRepeatabilityForBC(periodic_z_bcs, "z periodic"));
    PetscCall(CheckResidualRepeatabilityForBC(periodic_xy_bcs, "mixed x-y periodic"));
 
    PetscCall(PicurvMakeTempDir(profile_dir, sizeof(profile_dir)));
    PetscCall(PetscSNPrintf(profile_path, sizeof(profile_path), "%s/inlet.picslice", profile_dir));
    PetscCall(WriteNewtonPicSlice(profile_path));
    PetscCall(PetscSNPrintf(file_bcs, sizeof(file_bcs),
        "-Xi WALL noslip\n+Xi WALL noslip\n-Eta WALL noslip\n+Eta WALL noslip\n"
        "-Zeta INLET prescribed_velocity source_id=2 source_file=%s\n+Zeta OUTLET conservation\n", profile_path));
    PetscCall(CheckResidualRepeatabilityForBC(file_bcs, "file inlet/conservation outlet"));
    PetscCall(PicurvRemoveTempDir(profile_dir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies fixed, periodic-duplicate, and interior residual rows. */
static PetscErrorCode TestConstraintRows(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL, f = NULL;
    Cmpnts ***xa = NULL, ***fa = NULL, ***conditioned = NULL, ***rhs = NULL;
    MomentumNewtonKrylovContext ctx;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(geometric_periodic_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &f));
    PetscCall(VecSet(x, 0.25));
    PetscCall(DMDAVecGetArray(user->fda, x, &xa));
    if (user->info.xs == 0 && 2 >= user->info.ys && 2 < user->info.ys + user->info.ym &&
        2 >= user->info.zs && 2 < user->info.zs + user->info.zm) xa[2][2][0].x = 3.0;
    if (user->info.xs + user->info.xm == user->info.mx &&
        2 >= user->info.ys && 2 < user->info.ys + user->info.ym &&
        2 >= user->info.zs && 2 < user->info.zs + user->info.zm) xa[2][2][user->info.mx - 2].x = 1.25;
    PetscCall(DMDAVecRestoreArray(user->fda, x, &xa));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f, &ctx));
 
    PetscCall(DMDAVecGetArrayRead(user->fda, x, &xa));
    PetscCall(DMDAVecGetArrayRead(user->fda, f, &fa));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &conditioned));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Rhs, &rhs));
    if (user->info.xs == 0 && 2 >= user->info.ys && 2 < user->info.ys + user->info.ym &&
        2 >= user->info.zs && 2 < user->info.zs + user->info.zm) {
        PetscCall(PicurvAssertRealNear(1.75, fa[2][2][0].x, 1.0e-12,
                                       "periodic duplicate row must be Xdup-Xrep"));
    }
    if (user->info.ys == 0 && 2 >= user->info.xs && 2 < user->info.xs + user->info.xm &&
        2 >= user->info.zs && 2 < user->info.zs + user->info.zm) {
        PetscCall(PicurvAssertRealNear(xa[2][0][2].y - conditioned[2][0][2].y,
                                       fa[2][0][2].y, 1.0e-12,
                                       "fixed wall row must be X minus conditioned boundary value"));
    }
    if (2 >= user->info.xs && 2 < user->info.xs + user->info.xm &&
        2 >= user->info.ys && 2 < user->info.ys + user->info.ym &&
        2 >= user->info.zs && 2 < user->info.zs + user->info.zm) {
        PetscCall(PicurvAssertRealNear(-rhs[2][2][2].z, fa[2][2][2].z, 1.0e-12,
                                       "unconstrained interior row must retain the physical residual"));
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Rhs, &rhs));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &conditioned));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, f, &fa));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, x, &xa));
 
    PetscCall(VecDestroy(&f));
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Proves unit derivatives for every nonperiodic stored-row category and face. */
static PetscErrorCode TestFixedConstraintDerivativesAllFaces(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL;
    const PetscInt size[3] = {7, 7, 7};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecSet(x, 0.2));
 
    for (PetscInt axis = 0; axis < 3; ++axis) {
        PetscInt coord[3] = {2, 2, 2}, ri, rj, rk;
        PetscInt tangent = (axis + 1) % 3;
        MomentumRowType row;
 
        coord[axis] = 0;
        row = ClassifyMomentumRow(user, coord[0], coord[1], coord[2], axis, &ri, &rj, &rk);
        PetscCall(PicurvAssertIntEqual(MOM_ROW_FIXED_CONDITIONED, row,
                                       "negative face normal row classification"));
        PetscCall(CheckStoredDerivative(user, x, coord[0], coord[1], coord[2], axis,
            coord[0], coord[1], coord[2], axis, 1.0, 1.0e-8,
            "negative face normal fixed derivative"));
 
        row = ClassifyMomentumRow(user, coord[0], coord[1], coord[2], tangent, &ri, &rj, &rk);
        PetscCall(PicurvAssertIntEqual(MOM_ROW_FIXED_HOMOGENEOUS, row,
                                       "negative face tangential row classification"));
        PetscCall(CheckStoredDerivative(user, x, coord[0], coord[1], coord[2], tangent,
            coord[0], coord[1], coord[2], tangent, 1.0, 1.0e-8,
            "negative face tangential homogeneous derivative"));
 
        coord[axis] = size[axis] - 2;
        row = ClassifyMomentumRow(user, coord[0], coord[1], coord[2], axis, &ri, &rj, &rk);
        PetscCall(PicurvAssertIntEqual(MOM_ROW_FIXED_CONDITIONED, row,
                                       "positive physical normal row classification"));
        PetscCall(CheckStoredDerivative(user, x, coord[0], coord[1], coord[2], axis,
            coord[0], coord[1], coord[2], axis, 1.0, 1.0e-8,
            "positive physical normal fixed derivative"));
        row = ClassifyMomentumRow(user, coord[0], coord[1], coord[2], tangent, &ri, &rj, &rk);
        PetscCall(PicurvAssertIntEqual(MOM_ROW_PHYSICAL, row,
                                       "positive physical tangential row classification"));
 
        coord[axis] = size[axis] - 1;
        for (PetscInt component = 0; component < 3; ++component) {
            row = ClassifyMomentumRow(user, coord[0], coord[1], coord[2], component, &ri, &rj, &rk);
            PetscCall(PicurvAssertIntEqual(MOM_ROW_FIXED_HOMOGENEOUS, row,
                                           "positive dummy row classification"));
            PetscCall(CheckStoredDerivative(user, x, coord[0], coord[1], coord[2], component,
                coord[0], coord[1], coord[2], component, 1.0, 1.0e-8,
                "positive dummy homogeneous derivative"));
        }
    }
 
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Proves admitted inlet and outlet face-normal rows have unit self derivatives. */
static PetscErrorCode TestInletOutletConstraintDerivatives(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(NULL, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecCopy(user->Ucont, x));
    PetscCall(VecShift(x, 0.1));
    PetscCall(CheckStoredDerivative(user, x, 2, 2, 0, 2, 2, 2, 0, 2,
                                     1.0, 1.0e-8, "constant inlet fixed derivative"));
    /* The constant-velocity inlet imposes a value independent of X, so its
     * conditioned row F = X - cv has an exact unit self derivative.
     *
     * The conservation outlet is different: cv is the corrected outlet flux,
     * which the deterministic residual now reconstructs from the current X
     * (Ucat is seeded from X before the first outlet pass). Perturbing the
     * outlet-normal DOF therefore changes cv, so the self derivative is
     * 1 - dcv/dX and is strictly less than one. A self derivative of exactly
     * 1.0 here was an artifact of the pre-fix residual reading a stale
     * Cartesian state, i.e. an outlet correction decoupled from X. Assert the
     * derivative is (a) deterministic across independent evaluations -- the
     * residual-purity property -- and (b) reflects real conservation coupling
     * (0 < d < 1), rather than asserting a fixture-specific magic number. */
    {
        PetscReal d0 = 0.0, d1 = 0.0;
        PetscCall(MeasureStoredDerivative(user, x, 2, 2, user->info.mz - 2, 2,
                                          2, 2, user->info.mz - 2, 2, &d0));
        PetscCall(MeasureStoredDerivative(user, x, 2, 2, user->info.mz - 2, 2,
                                          2, 2, user->info.mz - 2, 2, &d1));
        PetscCall(PicurvAssertRealNear(d0, d1, 1.0e-9,
            "conservation outlet self derivative must be deterministic"));
        PetscCheck(d0 > 1.0e-3 && d0 < 1.0 - 1.0e-3, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
            "conservation outlet self derivative must reflect X-coupling (0<d<1), got %g.",
            (double)d0);
    }
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Checks one periodic configuration's endpoint derivatives on every component.
 * @param bcs Boundary text selecting the periodic axis.
 * @param axis Periodic axis index.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode CheckSingleAxisPeriodicDerivatives(const char *bcs, PetscInt axis)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL;
    PetscInt size[3], dup[3] = {2, 2, 2}, rep[3] = {2, 2, 2}, unrelated[3] = {3, 3, 3};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    size[0] = user->info.mx; size[1] = user->info.my; size[2] = user->info.mz;
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecSet(x, 0.15));
    for (PetscInt side = 0; side < 2; ++side) {
        dup[axis] = side == 0 ? 0 : size[axis] - 1;
        rep[axis] = side == 0 ? size[axis] - 2 : 1;
        for (PetscInt component = 0; component < 3; ++component) {
            PetscCall(CheckStoredDerivative(user, x, dup[0], dup[1], dup[2], component,
                dup[0], dup[1], dup[2], component, 1.0, 1.0e-8,
                "periodic duplicate self derivative"));
            PetscCall(CheckStoredDerivative(user, x, dup[0], dup[1], dup[2], component,
                rep[0], rep[1], rep[2], component, -1.0, 1.0e-8,
                "periodic representative derivative"));
            PetscCall(CheckStoredDerivative(user, x, dup[0], dup[1], dup[2], component,
                unrelated[0], unrelated[1], unrelated[2], (component + 1) % 3, 0.0, 1.0e-8,
                "periodic constraint unrelated derivative"));
        }
    }
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Proves single-, double-, triple-, and mixed-boundary periodic equations. */
static PetscErrorCode TestPeriodicConstraintDerivativesAndIntersections(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL, f = NULL;
    MomentumNewtonKrylovContext ctx;
    PetscScalar xdup, synced, residual;
 
    PetscFunctionBeginUser;
    PetscCall(CheckSingleAxisPeriodicDerivatives(periodic_x_bcs, 0));
    PetscCall(CheckSingleAxisPeriodicDerivatives(periodic_y_bcs, 1));
    PetscCall(CheckSingleAxisPeriodicDerivatives(periodic_z_bcs, 2));
 
    PetscCall(BuildNewtonFixture(periodic_xy_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &f));
    PetscCall(VecSet(x, 0.0));
    PetscCall(PerturbStoredValue(user, x, 0, 0, 2, 0, 3.0));
    PetscCall(PerturbStoredValue(user, x, user->info.mx - 2, user->info.my - 2, 2, 0, 1.25));
    PetscCall(VecCopy(x, user->Ucont));
    { const FieldId fields[] = {FIELD_ID_UCONT}; PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields)); }
    PetscCall(GetStoredValue(user, x, 0, 0, 2, 0, &xdup));
    PetscCall(GetStoredValue(user, user->Ucont, 0, 0, 2, 0, &synced));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f, &ctx));
    PetscCall(GetStoredValue(user, f, 0, 0, 2, 0, &residual));
    PetscCall(PicurvAssertRealNear(PetscRealPart(xdup - synced), PetscRealPart(residual), 1.0e-12,
                                   "doubly periodic edge must use production synchronized representative"));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 2, 0, 0, 0, 2, 0, 1.0, 1.0e-8,
                                     "doubly periodic edge self derivative"));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 2, 0,
        user->info.mx - 2, user->info.my - 2, 2, 0, -1.0, 1.0e-8,
        "doubly periodic edge representative derivative"));
    PetscCall(VecDestroy(&f)); PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
 
    tmpdir[0] = '\0';
    PetscCall(BuildNewtonFixture(periodic_xyz_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x)); PetscCall(VecSet(x, 0.1));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 0, 2, 0, 0, 0, 2, 1.0, 1.0e-8,
                                     "fully periodic corner self derivative"));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 0, 2,
        user->info.mx - 2, user->info.my - 2, user->info.mz - 2, 2,
        -1.0, 1.0e-8, "fully periodic corner representative derivative"));
    PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
 
    tmpdir[0] = '\0';
    PetscCall(BuildNewtonFixture(periodic_x_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs)); PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecSet(x, 0.1));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 2, 1, 0, 0, 2, 1, 1.0, 1.0e-8,
                                     "periodic-wall intersection self derivative"));
    PetscCall(CheckStoredDerivative(user, x, 0, 0, 2, 1,
        user->info.mx - 2, 0, 2, 1, -1.0, 1.0e-8,
        "periodic-wall intersection representative derivative"));
    PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Compares PETSc's matrix-free action with direct differencing. */
static PetscErrorCode TestMatrixFreeDerivative(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    SNES snes = NULL;
    Mat J = NULL;
    Vec x = NULL, xp = NULL, f0 = NULL, fp = NULL, v = NULL, jv = NULL, fd = NULL;
    PetscReal h = 0.0, error = 0.0, scale = 0.0;
    MomentumNewtonKrylovContext ctx;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(NULL, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &xp));
    PetscCall(VecDuplicate(user->Ucont, &f0));
    PetscCall(VecDuplicate(user->Ucont, &fp));
    PetscCall(VecDuplicate(user->Ucont, &v));
    PetscCall(VecDuplicate(user->Ucont, &jv));
    PetscCall(VecDuplicate(user->Ucont, &fd));
    PetscCall(VecCopy(user->Ucont, x));
    PetscCall(VecShift(x, 0.05));
    PetscCall(VecSet(v, 0.5));
    ctx.user = user;
 
    PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    PetscCall(SNESSetDM(snes, user->fda));
    PetscCall(SNESSetFunction(snes, f0, MomentumNewtonKrylov_FormResidual, &ctx));
    PetscCall(MatCreateSNESMF(snes, &J));
    PetscCall(MomentumNewtonKrylov_FormResidual(snes, x, f0, &ctx));
    PetscCall(MatMFFDSetBase(J, x, f0));
    PetscCall(MatMult(J, v, jv));
    PetscCall(MatMFFDGetH(J, &h));
    PetscCall(VecWAXPY(xp, h, v, x));
    PetscCall(MomentumNewtonKrylov_FormResidual(snes, xp, fp, &ctx));
    PetscCall(VecWAXPY(fd, -1.0, f0, fp));
    PetscCall(VecScale(fd, 1.0 / h));
    PetscCall(VecAXPY(fd, -1.0, jv));
    PetscCall(VecNorm(fd, NORM_2, &error));
    PetscCall(VecNorm(jv, NORM_2, &scale));
    PetscCall(PicurvAssertBool((PetscBool)(error <= 1.0e-9 * PetscMax(1.0, scale)),
                               "matrix-free Jv must match direct differencing"));
 
    PetscCall(VecDestroy(&fd));
    PetscCall(VecDestroy(&jv));
    PetscCall(VecDestroy(&v));
    PetscCall(VecDestroy(&fp));
    PetscCall(VecDestroy(&f0));
    PetscCall(VecDestroy(&xp));
    PetscCall(VecDestroy(&x));
    PetscCall(MatDestroy(&J));
    PetscCall(SNESDestroy(&snes));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Builds a compact all-wall operator fixture through real boundary handlers.
 * @param simCtx Returned simulation context.
 * @param user Returned block context.
 * @param x_periodic Whether the x faces use geometric periodicity.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode BuildMinimalWallOperatorFixture(SimCtx **simCtx, UserCtx **user,
                                                       PetscBool x_periodic)
{
    PetscFunctionBeginUser;
    /* Request the production-width (3) DMDA stencil so the complete RHS is safe
       across MPI partitions; boundary metadata below still selects physical walls. */
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(
        simCtx, user, 6, 6, 6, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    (*simCtx)->i_periodic = x_periodic ? 1 : 0;
    (*simCtx)->j_periodic = (*simCtx)->k_periodic = 0;
    (*simCtx)->mom_solver_type = MOMENTUM_SOLVER_NEWTON_KRYLOV;
    (*simCtx)->invicid = 1;
    (*simCtx)->dt = 0.1;
    (*simCtx)->step = 1;
    (*simCtx)->StartStep = 0;
    PetscCall(VecSet((*user)->Ucont, 0.0));
    PetscCall(VecSet((*user)->Ucont_o, 0.0));
    PetscCall(VecSet((*user)->Ucont_rm1, 0.0));
    for (PetscInt face = 0; face < 6; ++face) {
        PetscBool periodic_face = (PetscBool)(x_periodic &&
            (face == BC_FACE_NEG_X || face == BC_FACE_POS_X));
        (*user)->boundary_faces[face].face_id = (BCFace)face;
        (*user)->boundary_faces[face].mathematical_type = periodic_face ? PERIODIC : WALL;
        (*user)->boundary_faces[face].handler_type = periodic_face ?
            BC_HANDLER_PERIODIC_GEOMETRIC : BC_HANDLER_WALL_NOSLIP;
        PetscCall(BoundaryCondition_Create((*user)->boundary_faces[face].handler_type,
                                            &(*user)->boundary_faces[face].handler));
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief A Newton solve with free-slip walls converges, preconditioned and not.
 *
 * Free-slip walls are supported by Newton Krylov from the start: the matrix-free action
 * carries the reflected boundary through the residual, and the point-block
 * preconditioner approximates the operator without wall-specific terms.
 */
static PetscErrorCode TestFreeSlipWallsSolve(void)
{
    const char *models[] = {"none", "frozen_momentum_jacobian"};
    const char *structures[] = {"none", "point_block"};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(models) / sizeof(models[0]); ++n) {
        SimCtx *simCtx = NULL;
        UserCtx *user = NULL;
        char tmpdir[PETSC_MAX_PATH_LEN] = "";
 
        PetscCall(BuildNewtonFixture(free_slip_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-4"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-6"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model", models[n]));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", structures[n]));
        PetscCall(VecSet(user->Ucont, 0.1));
        PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
        PetscCall(PicurvAssertBool(simCtx->mom_last_converged, "a solve with free-slip walls must converge"));
        PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
        PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_model"));
        PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_structure"));
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief A TwoD solve converges, preconditioned and not, and holds the masked component.
 *
 * ComputeRHS() zeroes the masked component's spatial terms, so its residual row is the
 * time derivative alone. With the same value at every history level, BDF2's solution of
 * that row is the value itself, so a nonzero xi flux must survive the solve unchanged, as
 * the Picard solver's skipped row keeps it.
 */
static PetscErrorCode TestTwoDSolveHoldsMaskedComponent(void)
{
    const char *models[] = {"none", "frozen_momentum_jacobian"};
    const char *structures[] = {"none", "point_block"};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(models) / sizeof(models[0]); ++n) {
        SimCtx *simCtx = NULL;
        UserCtx *user = NULL;
        Cmpnts ***ucont = NULL;
        PetscReal drift = 0.0, global_drift = 0.0;
        char tmpdir[PETSC_MAX_PATH_LEN] = "";
 
        PetscCall(BuildNewtonFixture(periodic_x_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
        simCtx->TwoD = 1;
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-8"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-10"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model", models[n]));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", structures[n]));
        PetscCall(VecSet(user->Ucont, 0.0));
        PetscCall(VecStrideSet(user->Ucont, 0, 0.1));
        PetscCall(VecStrideSet(user->Ucont, 2, 0.05));
        PetscCall(VecCopy(user->Ucont, user->Ucont_o));
        PetscCall(VecCopy(user->Ucont, user->Ucont_rm1));
        PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
        PetscCall(PicurvAssertBool(simCtx->mom_last_converged, "a TwoD solve must converge"));
 
        PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
        for (PetscInt k = PetscMax(user->info.zs, 1); k < PetscMin(user->info.zs + user->info.zm, user->info.mz - 2); ++k)
            for (PetscInt j = PetscMax(user->info.ys, 1); j < PetscMin(user->info.ys + user->info.ym, user->info.my - 2); ++j)
                for (PetscInt i = PetscMax(user->info.xs, 1); i < PetscMin(user->info.xs + user->info.xm, user->info.mx - 1); ++i)
                    drift = PetscMax(drift, PetscAbsReal(ucont[k][j][i].x - 0.1));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
        PetscCallMPI(MPI_Allreduce(&drift, &global_drift, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
        PetscCall(PicurvAssertRealNear(0.0, global_drift, 1.0e-12,
                                       "the masked component must keep its time-constant value"));
        PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
        PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_model"));
        PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_structure"));
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Forms the complete direct FD Jacobian, checks every row, and compares MFFD actions.
 */
static PetscErrorCode TestWholeOperatorDirectJacobian(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    SNES snes = NULL;
    Mat J = NULL;
    Vec x = NULL, xp = NULL, f0 = NULL, fp = NULL, column = NULL, square = NULL;
    Vec row_norm_sq = NULL, v[2] = {NULL, NULL}, dense_v[2] = {NULL, NULL};
    Vec mffd_v = NULL, error_vec = NULL;
    PetscInt n_global, lo, hi;
    PetscReal min_row_sq = 0.0, error = 0.0, reference = 0.0;
    const PetscReal epsilon = 1.0e-7;
    MomentumNewtonKrylovContext ctx;
 
    PetscFunctionBeginUser;
    PetscCall(BuildMinimalWallOperatorFixture(&simCtx, &user, PETSC_FALSE));
    PetscCall(VecDuplicate(user->Ucont, &x)); PetscCall(VecSet(x, 0.0));
    PetscCall(VecDuplicate(x, &xp)); PetscCall(VecDuplicate(x, &f0));
    PetscCall(VecDuplicate(x, &fp)); PetscCall(VecDuplicate(x, &column));
    PetscCall(VecDuplicate(x, &square)); PetscCall(VecDuplicate(x, &row_norm_sq));
    PetscCall(VecDuplicate(x, &v[0])); PetscCall(VecDuplicate(x, &v[1]));
    PetscCall(VecDuplicate(x, &dense_v[0])); PetscCall(VecDuplicate(x, &dense_v[1]));
    PetscCall(VecDuplicate(x, &mffd_v)); PetscCall(VecDuplicate(x, &error_vec));
    PetscCall(VecZeroEntries(row_norm_sq)); PetscCall(VecZeroEntries(dense_v[0]));
    PetscCall(VecZeroEntries(dense_v[1]));
    PetscCall(VecGetSize(x, &n_global));
    PetscCall(VecGetOwnershipRange(x, &lo, &hi));
    for (PetscInt which = 0; which < 2; ++which) {
        PetscScalar *a = NULL;
        PetscCall(VecGetArray(v[which], &a));
        for (PetscInt local = 0; local < hi - lo; ++local) {
            PetscInt global = lo + local;
            a[local] = which == 0 ? (PetscScalar)(1.0 + 0.05 * (global % 9))
                                  : (PetscScalar)(((global % 2) ? -1.0 : 1.0) * (0.5 + 0.03 * (global % 7)));
        }
        PetscCall(VecRestoreArray(v[which], &a));
    }
 
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    for (PetscInt col = 0; col < n_global; ++col) {
        PetscScalar coeff[2] = {
            (PetscScalar)(1.0 + 0.05 * (col % 9)),
            (PetscScalar)(((col % 2) ? -1.0 : 1.0) * (0.5 + 0.03 * (col % 7)))
        };
        PetscCall(VecCopy(x, xp));
        if (col >= lo && col < hi) PetscCall(VecSetValue(xp, col, epsilon, ADD_VALUES));
        PetscCall(VecAssemblyBegin(xp)); PetscCall(VecAssemblyEnd(xp));
        PetscCall(MomentumNewtonKrylov_FormResidual(NULL, xp, fp, &ctx));
        PetscCall(VecWAXPY(column, -1.0, f0, fp));
        PetscCall(VecScale(column, 1.0 / epsilon));
        PetscCall(VecPointwiseMult(square, column, column));
        PetscCall(VecAXPY(row_norm_sq, 1.0, square));
        PetscCall(VecAXPY(dense_v[0], coeff[0], column));
        PetscCall(VecAXPY(dense_v[1], coeff[1], column));
    }
    PetscCall(VecMin(row_norm_sq, NULL, &min_row_sq));
    PetscCheck(min_row_sq > 0.5, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "Complete Newton Jacobian contains an unexplained zero/weak row (min squared norm=%g).",
               (double)min_row_sq);
    PetscCall(CheckStoredDerivative(user, x, 2, 2, 2, 0, 2, 2, 2, 0,
                                     10.0, 1.0e-5, "interior BDF1 temporal diagonal"));
 
    PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    PetscCall(SNESSetDM(snes, user->fda));
    PetscCall(SNESSetFunction(snes, f0, MomentumNewtonKrylov_FormResidual, &ctx));
    PetscCall(MatCreateSNESMF(snes, &J));
    PetscCall(MatMFFDSetBase(J, x, f0));
    for (PetscInt which = 0; which < 2; ++which) {
        PetscCall(MatMult(J, v[which], mffd_v));
        PetscCall(VecWAXPY(error_vec, -1.0, dense_v[which], mffd_v));
        PetscCall(VecNorm(error_vec, NORM_2, &error));
        PetscCall(VecNorm(dense_v[which], NORM_2, &reference));
        PetscCheck(error <= 2.0e-5 * PetscMax(1.0, reference), PETSC_COMM_WORLD, PETSC_ERR_PLIB,
                   "Independent dense FD Jv differs from PETSc MFFD action %" PetscInt_FMT ": error=%g reference=%g.",
                   which, (double)error, (double)reference);
    }
 
    for (PetscInt which = 0; which < 2; ++which) {
        const PetscReal steps[3] = {1.0e-4, 1.0e-6, 1.0e-8};
        PetscReal best = PETSC_MAX_REAL;
        for (PetscInt s = 0; s < 3; ++s) {
            PetscCall(VecWAXPY(xp, steps[s], v[which], x));
            PetscCall(MomentumNewtonKrylov_FormResidual(NULL, xp, fp, &ctx));
            PetscCall(VecWAXPY(column, -1.0, f0, fp));
            PetscCall(VecScale(column, 1.0 / steps[s]));
            PetscCall(VecAXPY(column, -1.0, dense_v[which]));
            PetscCall(VecNorm(column, NORM_2, &error));
            best = PetscMin(best, error);
        }
        PetscCall(VecNorm(dense_v[which], NORM_2, &reference));
        PetscCheck(best <= 2.0e-5 * PetscMax(1.0, reference), PETSC_COMM_WORLD, PETSC_ERR_PLIB,
                   "Direct directional differences show no accuracy plateau for vector %" PetscInt_FMT " (best=%g).",
                   which, (double)best);
    }
 
    PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
    PetscCall(VecDestroy(&error_vec)); PetscCall(VecDestroy(&mffd_v));
    PetscCall(VecDestroy(&dense_v[1])); PetscCall(VecDestroy(&dense_v[0]));
    PetscCall(VecDestroy(&v[1])); PetscCall(VecDestroy(&v[0]));
    PetscCall(VecDestroy(&row_norm_sq)); PetscCall(VecDestroy(&square));
    PetscCall(VecDestroy(&column)); PetscCall(VecDestroy(&fp)); PetscCall(VecDestroy(&f0));
    PetscCall(VecDestroy(&xp)); PetscCall(VecDestroy(&x));
    PetscCall(BoundarySystem_Destroy(user));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Audits every row of a complete operator containing periodic duplicates. */
static PetscErrorCode TestPeriodicOperatorHasNoZeroRows(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    Vec x = NULL, xp = NULL, f0 = NULL, fp = NULL, column = NULL;
    Vec square = NULL, row_norm_sq = NULL;
    PetscInt n_global, lo, hi;
    PetscReal min_row_sq = 0.0;
    const PetscReal epsilon = 1.0e-7;
    MomentumNewtonKrylovContext ctx;
 
    PetscFunctionBeginUser;
    PetscCall(BuildMinimalWallOperatorFixture(&simCtx, &user, PETSC_TRUE));
    PetscCall(VecDuplicate(user->Ucont, &x)); PetscCall(VecSet(x, 0.0));
    PetscCall(VecDuplicate(x, &xp)); PetscCall(VecDuplicate(x, &f0));
    PetscCall(VecDuplicate(x, &fp)); PetscCall(VecDuplicate(x, &column));
    PetscCall(VecDuplicate(x, &square)); PetscCall(VecDuplicate(x, &row_norm_sq));
    PetscCall(VecZeroEntries(row_norm_sq));
    PetscCall(VecGetSize(x, &n_global)); PetscCall(VecGetOwnershipRange(x, &lo, &hi));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f0, &ctx));
    for (PetscInt col = 0; col < n_global; ++col) {
        PetscCall(VecCopy(x, xp));
        if (col >= lo && col < hi) PetscCall(VecSetValue(xp, col, epsilon, ADD_VALUES));
        PetscCall(VecAssemblyBegin(xp)); PetscCall(VecAssemblyEnd(xp));
        PetscCall(MomentumNewtonKrylov_FormResidual(NULL, xp, fp, &ctx));
        PetscCall(VecWAXPY(column, -1.0, f0, fp)); PetscCall(VecScale(column, 1.0 / epsilon));
        PetscCall(VecPointwiseMult(square, column, column));
        PetscCall(VecAXPY(row_norm_sq, 1.0, square));
    }
    PetscCall(VecMin(row_norm_sq, NULL, &min_row_sq));
    PetscCheck(min_row_sq > 0.5, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "Periodic Newton Jacobian contains a zero/weak row (min squared norm=%g).",
               (double)min_row_sq);
    PetscCall(VecDestroy(&row_norm_sq)); PetscCall(VecDestroy(&square));
    PetscCall(VecDestroy(&column)); PetscCall(VecDestroy(&fp)); PetscCall(VecDestroy(&f0));
    PetscCall(VecDestroy(&xp)); PetscCall(VecDestroy(&x));
    PetscCall(BoundarySystem_Destroy(user));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Exercises a converged solve, a committed unconverged solve, a non-finite rollback, and cleanup. */
static PetscErrorCode TestSmallSolveAndRollback(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec entry = NULL, delta = NULL;
    PetscErrorCode solve_ierr;
    PetscReal norm = 0.0;
    const FieldId fields[] = {FIELD_ID_UCONT};
    char summary_path[PETSC_MAX_PATH_LEN];
    char history_path[PETSC_MAX_PATH_LEN];
    char linear_history_path[PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-4"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-6"));
    simCtx->mom_nk_monitor_history = PETSC_TRUE;
    PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
    PetscCall(PicurvAssertBool(simCtx->mom_last_converged, "small Newton solve must converge"));
    PetscCall(PicurvAssertBool((PetscBool)(user->Rhs != NULL && user->momentumNewtonKrylov != NULL),
                               "successful solve must keep its workspace for the next step"));
    PetscCall(PetscSNPrintf(summary_path, sizeof(summary_path),
                            "%s/Momentum_Solver_Newton_Krylov_Summary_Block_0.csv", simCtx->log_dir));
    PetscCall(PetscSNPrintf(history_path, sizeof(history_path),
                            "%s/Momentum_Solver_Newton_Krylov_History_Block_0.csv", simCtx->log_dir));
    PetscCall(PetscSNPrintf(linear_history_path, sizeof(linear_history_path),
                            "%s/Momentum_Solver_Newton_Krylov_Linear_History_Block_0.csv",
                            simCtx->log_dir));
    PetscCall(AssertNewtonLog(summary_path, 1, ",Newton Krylov,", ",committed\n"));
    PetscCall(AssertNewtonLog(history_path, -2, "step,block,newton,nonlinear_norm", ",0,0,"));
    PetscCall(AssertNewtonLog(linear_history_path, -2, "krylov,requested_rtol,reported_residual_norm", ",0,0,0,"));
 
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    tmpdir[0] = '\0';
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecSet(user->Ucont, 0.2));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(VecDuplicate(user->Ucont, &entry));
    PetscCall(VecDuplicate(user->Ucont, &delta));
    PetscCall(VecCopy(user->Ucont, entry));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "0"));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PicurvAssertIntEqual(0, solve_ierr,
                                   "a finite residual no larger than at entry must commit the step"));
    PetscCall(PicurvAssertBool((PetscBool)!simCtx->mom_last_converged,
                               "a committed SNES failure must be reported as not converged"));
    PetscCall(VecWAXPY(delta, -1.0, entry, user->Ucont));
    PetscCall(VecNorm(delta, NORM_INFINITY, &norm));
    PetscCall(PicurvAssertRealNear(0.0, norm, 1.0e-12,
                                   "a zero-iteration failure must commit the unchanged entry state"));
    PetscCall(PicurvAssertBool((PetscBool)(user->Rhs != NULL && user->momentumNewtonKrylov != NULL),
                               "committed unconverged solve must keep its workspace for the next step"));
    PetscCall(PetscSNPrintf(summary_path, sizeof(summary_path),
                            "%s/Momentum_Solver_Newton_Krylov_Summary_Block_0.csv", simCtx->log_dir));
    PetscCall(AssertNewtonLog(summary_path, 1, ",DIVERGED_", ",committed_unconverged\n"));
 
    /* A non-finite residual is the one failure that must still roll back and stop the run.
       Options are read when the solve context is built, so release it to apply the new limit. */
    PetscCall(MomentumSolver_NewtonKrylov_Destroy(user));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
    PetscCall(VecSet(user->Ucont, PETSC_MAX_REAL));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(VecCopy(user->Ucont, entry));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_CONV_FAILED, solve_ierr,
                                   "a non-finite residual must report PETSC_ERR_CONV_FAILED"));
    PetscCall(PicurvAssertBool((PetscBool)!simCtx->mom_last_converged,
                               "a rolled-back solve must be reported as not converged"));
    PetscCall(VecWAXPY(delta, -1.0, entry, user->Ucont));
    PetscCall(VecNorm(delta, NORM_INFINITY, &norm));
    PetscCall(PicurvAssertRealNear(0.0, norm, 1.0e-12,
                                   "failed Newton solve must restore the canonical entry state"));
    PetscCall(PicurvAssertBool((PetscBool)(user->Rhs == NULL), "failed solve must release Rhs"));
    PetscCall(AssertNewtonLog(summary_path, 2, ",DIVERGED_", ",rolled_back\n"));
 
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_rtol"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_max_it"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_ksp_rtol"));
    PetscCall(VecDestroy(&delta));
    PetscCall(VecDestroy(&entry));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies the six-wall zero-velocity case logs zero Newton/Krylov work. */
static PetscErrorCode TestZeroIterationStructuredLogging(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    char summary_path[PETSC_MAX_PATH_LEN];
    const FieldId fields[] = {FIELD_ID_UCONT};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecZeroEntries(user->Ucont));
    PetscCall(VecZeroEntries(user->Ucont_o));
    PetscCall(VecZeroEntries(user->Ucont_rm1));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
    PetscCall(PetscSNPrintf(summary_path, sizeof(summary_path),
                            "%s/Momentum_Solver_Newton_Krylov_Summary_Block_0.csv", simCtx->log_dir));
    PetscCall(AssertNewtonLog(summary_path, 1, ",0,1,0,",
                             ",0.0000000000000000e+00,committed\n"));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Exercises the straight-duct BDF1 startup path used by flat_channel.
 *
 * The conservation outlet consumes lUcat during its first boundary pass.  This
 * test deliberately evaluates the callback at the initialized state before
 * installing SNES, then completes the first nonlinear solve with each shipped
 * preconditioner.  It catches a missing Ucont -> Ucat -> lUcat seed as a
 * non-finite initial residual rather than hiding it behind later MFFD work.
 */
static PetscErrorCode CheckFlatChannelStartup(PetscBool use_point_block)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec x = NULL, f = NULL;
    PetscReal initial_norm = 0.0;
    MomentumNewtonKrylovContext ctx = {0};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(NULL, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &f));
    PetscCall(VecCopy(user->Ucont, x));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f, &ctx));
    PetscCall(VecNorm(f, NORM_2, &initial_norm));
    PetscCheck(!PetscIsInfOrNanReal(initial_norm), PETSC_COMM_WORLD, PETSC_ERR_FP,
               "flat-channel BDF1 initial residual is non-finite (%g) with %s.",
               (double)initial_norm,
               use_point_block ? "frozen-momentum point-block" : "PCNONE");
    PetscCall(VecDestroy(&f));
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-4"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-6"));
    if (use_point_block) {
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model",
                                       "frozen_momentum_jacobian"));
        PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure",
                                       "point_block"));
    }
    PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
    PetscCall(PicurvAssertBool(simCtx->mom_last_converged,
                               "flat-channel BDF1 Newton solve must converge"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_model"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_structure"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_rtol"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_max_it"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_ksp_rtol"));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Guards flat_channel's initial BDF1 residual and both shipped NK PCs. */
static PetscErrorCode TestFlatChannelStartup(void)
{
    PetscFunctionBeginUser;
    PetscCall(CheckFlatChannelStartup(PETSC_FALSE));
    PetscCall(CheckFlatChannelStartup(PETSC_TRUE));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies restarted Newton solves with both supported preconditioners. */
static PetscErrorCode TestRestartAndContinuationSolve(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    const FieldId fields[] = {FIELD_ID_UCONT};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecSet(user->Ucont, 0.2));
    PetscCall(VecZeroEntries(user->Ucont_o));
    PetscCall(VecZeroEntries(user->Ucont_rm1));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-4"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "20"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-6"));
 
    /* AdvanceSimulation() solves StartStep+1 first.  Use a nonzero checkpoint
     * state so that SNES takes a Newton/Krylov path rather than accepting a
     * trivial residual. */
    simCtx->StartStep = 7;
    simCtx->step = simCtx->StartStep + 1;
    PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
    PetscCall(PicurvAssertBool(simCtx->mom_last_converged,
                               "PCNONE Newton Krylov checkpoint restart must converge"));
 
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    tmpdir[0] = '\0';
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecSet(user->Ucont, 0.2));
    PetscCall(VecZeroEntries(user->Ucont_o));
    PetscCall(VecZeroEntries(user->Ucont_rm1));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model",
                                   "frozen_momentum_jacobian"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "point_block"));
 
    simCtx->continueMode = PETSC_TRUE;
    simCtx->StartStep = 8;
    simCtx->step = simCtx->StartStep + 1;
    PetscCall(MomentumSolver_NewtonKrylov(user, NULL, NULL));
    PetscCall(PicurvAssertBool(simCtx->mom_last_converged,
                               "frozen-momentum Newton Krylov --continue solve must converge"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_model"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_preconditioner_structure"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_rtol"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_snes_max_it"));
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_ksp_rtol"));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Confirms unsupported features fail before workspace allocation. */
static PetscErrorCode TestUnsupportedConfigurationFailsBeforeAllocation(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    PetscErrorCode solve_ierr;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    {
        /* The gradient (Clark) model and wall functions are no longer rejected: both
           reach the solver through the shared residual, which this solver already
           evaluates. Neither is verified under Newton-Krylov yet, so neither is
           claimed as supported; they are simply no longer refused. */
        PetscInt *unsupported_flags[] = {
            &simCtx->immersed, &simCtx->movefsi, &simCtx->rotatefsi,
            &simCtx->moveframe, &simCtx->rotateframe
        };
        for (size_t flag = 0; flag < sizeof(unsupported_flags) / sizeof(unsupported_flags[0]); ++flag) {
            *unsupported_flags[flag] = 1;
            PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
            solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
            PetscCall(PetscPopErrorHandler());
            PetscCall(PicurvAssertBool((PetscBool)(solve_ierr != PETSC_SUCCESS),
                                       "unsupported Newton feature flag must fail"));
            PetscCall(PicurvAssertBool((PetscBool)(user->Rhs == NULL),
                                       "feature validation must precede workspace allocation"));
            *unsupported_flags[flag] = 0;
        }
    }
    simCtx->block_number = 2;
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(solve_ierr != PETSC_SUCCESS), "multiblock must fail"));
    simCtx->block_number = 1;
    /* Solid cells are no longer a rejection. The residual already zeroes their rows, so
       they carry no unknown; the solver now gives them the same identity treatment as
       any other row without one, instead of refusing the whole solve. A fully masked
       domain is the degenerate case of that: every row is an identity row and the
       system is trivially consistent. This checks the structural contract only; the
       masked path has no physics verification because nothing in the tree currently
       populates Nvert with a nonzero value. */
    PetscCall(VecSet(user->Nvert, 1.0));
    PetscCall(VecSet(user->lNvert, 1.0));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PicurvAssertBool((PetscBool)(solve_ierr == PETSC_SUCCESS),
                               "a fully masked domain must resolve to identity rows, not a rejection"));
    /* That solve kept a context built for the all-solid mask; release it so the remaining
       checks see the fixture's unallocated workspace again. */
    PetscCall(MomentumSolver_NewtonKrylov_Destroy(user));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(VecSet(user->Nvert, 0.0));
    PetscCall(VecSet(user->lNvert, 0.0));
    PetscCall(VecSet(user->Nvert, 0.0));
    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;
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertBool((PetscBool)(solve_ierr != PETSC_SUCCESS),
                               "driven constant-flux controller must fail"));
    PetscCall(PicurvAssertBool((PetscBool)(user->Rhs == NULL),
                               "all validation failures must precede workspace allocation"));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies cleanup and rollback after an options failure following asset creation. */
static PetscErrorCode TestPostAllocationFailureCleanup(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    Vec entry = NULL, delta = NULL;
    PetscErrorCode solve_ierr;
    PetscReal norm = 0.0;
    const FieldId fields[] = {FIELD_ID_UCONT};
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecSet(user->Ucont, 0.2));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(VecDuplicate(user->Ucont, &entry)); PetscCall(VecCopy(user->Ucont, entry));
    PetscCall(VecDuplicate(user->Ucont, &delta));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_pc_type", "jacobi"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    solve_ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PetscOptionsClearValue(NULL, "-mom_nk_pc_type"));
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_SUP, solve_ierr,
                                   "non-PCNONE option must fail after setup"));
    PetscCall(PicurvAssertBool((PetscBool)(user->Rhs == NULL),
                               "post-allocation failure must destroy Rhs"));
    PetscCall(VecWAXPY(delta, -1.0, entry, user->Ucont));
    PetscCall(VecNorm(delta, NORM_INFINITY, &norm));
    PetscCall(PicurvAssertRealNear(0.0, norm, 1.0e-12,
                                   "post-allocation failure must restore canonical entry"));
    PetscCall(VecDestroy(&delta)); PetscCall(VecDestroy(&entry));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies finalized application-owned linearization option parsing. */
static PetscErrorCode TestLinearizationConfigParsing(void)
{
    MomentumNewtonJacobian jacobian = {0};
    MomentumPreconditionerDescription description = {0};
    PetscErrorCode config_ierr;
    const char *option_names[] = {
        "-mom_nk_jacobian_type",
        "-mom_nk_jacobian_fd_mode",
        "-mom_nk_preconditioner_model",
        "-mom_nk_preconditioner_structure"
    };
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(option_names) / sizeof(option_names[0]); ++n)
        PetscCall(PetscOptionsClearValue(NULL, option_names[n]));
 
    PetscCall(MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description));
    PetscCall(PicurvAssertIntEqual(MOM_NK_JACOBIAN_FINITE_DIFFERENCE, jacobian.type,
                                   "default Jacobian type must be finite difference"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_FD_MODE_MATRIX_FREE,
                                   jacobian.finite_difference_mode,
                                   "default finite-difference mode must be matrix free"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_MODEL_NONE, description.model,
                                   "default preconditioner model must be none"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_STRUCTURE_NONE, description.structure,
                                   "default preconditioner structure must be none"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_type", "finite_difference"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_fd_mode", "matrix_free"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model", "none"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "none"));
    PetscCall(MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description));
    PetscCall(PicurvAssertIntEqual(MOM_NK_JACOBIAN_FINITE_DIFFERENCE, jacobian.type,
                                   "explicit baseline Jacobian type must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_FD_MODE_MATRIX_FREE,
                                   jacobian.finite_difference_mode,
                                   "explicit baseline finite-difference mode must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_MODEL_NONE, description.model,
                                   "explicit baseline preconditioner model must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_STRUCTURE_NONE, description.structure,
                                   "explicit baseline preconditioner structure must parse"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model",
                                  "frozen_momentum_jacobian"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "point_block"));
    PetscCall(MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description));
    PetscCall(PicurvAssertIntEqual(MOM_NK_JACOBIAN_FINITE_DIFFERENCE, jacobian.type,
                                   "explicit Jacobian type must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_FD_MODE_MATRIX_FREE,
                                   jacobian.finite_difference_mode,
                                   "explicit finite-difference mode must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN,
                                   description.model, "frozen model must parse"));
    PetscCall(PicurvAssertIntEqual(MOM_NK_PC_STRUCTURE_POINT_BLOCK,
                                   description.structure, "point-block structure must parse"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_type", "analytic"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONG, config_ierr,
                                   "unsupported Jacobian type must fail"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_type", "finite_difference"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_fd_mode", "colored_sparse"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONG, config_ierr,
                                   "unsupported finite-difference mode must fail"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_jacobian_fd_mode", "matrix_free"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "none"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_SUP, config_ierr,
                                   "frozen model without point block must fail"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model", "none"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "point_block"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_SUP, config_ierr,
                                   "none model with point block must fail"));
 
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model", "diagonal"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "none"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONG, config_ierr,
                                   "unsupported preconditioner model must fail"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_model",
                                  "frozen_momentum_jacobian"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_preconditioner_structure", "line"));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    config_ierr = MomentumNewtonKrylov_ReadLinearizationConfig(&jacobian, &description);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONG, config_ierr,
                                   "unsupported preconditioner structure must fail"));
 
    for (size_t n = 0; n < sizeof(option_names) / sizeof(option_names[0]); ++n)
        PetscCall(PetscOptionsClearValue(NULL, option_names[n]));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Converts an in-domain or periodic-ghost DMDA stencil to PETSc ordering. */
static PetscErrorCode TestStencilToGlobal(UserCtx *user, MatStencil stencil,
                                          PetscInt *global_index)
{
    ISLocalToGlobalMapping local_to_global = NULL;
    PetscInt ghost_starts[3], ghost_sizes[3], local_index;
 
    PetscFunctionBeginUser;
    PetscCall(DMDAGetGhostCorners(user->fda,
                                  &ghost_starts[0], &ghost_starts[1], &ghost_starts[2],
                                  &ghost_sizes[0], &ghost_sizes[1], &ghost_sizes[2]));
    PetscCall(PicurvAssertBool((PetscBool)(
        stencil.i >= ghost_starts[0] && stencil.i < ghost_starts[0] + ghost_sizes[0] &&
        stencil.j >= ghost_starts[1] && stencil.j < ghost_starts[1] + ghost_sizes[1] &&
        stencil.k >= ghost_starts[2] && stencil.k < ghost_starts[2] + ghost_sizes[2]),
        "test stencil must lie in the local DMDA ghost region"));
    local_index = stencil.c + 3 * (
        (stencil.i - ghost_starts[0]) + ghost_sizes[0] * (
        (stencil.j - ghost_starts[1]) + ghost_sizes[1] *
        (stencil.k - ghost_starts[2])));
    PetscCall(DMGetLocalToGlobalMapping(user->fda, &local_to_global));
    PetscCall(ISLocalToGlobalMappingApply(local_to_global, 1, &local_index,
                                          global_index));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Reads one DMDA-stencil matrix entry through collective basis vectors. */
static PetscErrorCode PreconditionerMatrixStencilEntry(UserCtx *user,
                                       Mat preconditioning_matrix, MatStencil row,
                                       MatStencil col, PetscScalar *value)
{
    Vec column_basis = NULL, row_basis = NULL, product = NULL;
    AO ao = NULL;
    PetscInt mx, my, row_index, col_index;
    PetscFunctionBeginUser;
    PetscCall(DMDAGetInfo(user->fda, NULL, &mx, &my, NULL, NULL, NULL, NULL,
                          NULL, NULL, NULL, NULL, NULL, NULL));
    row_index = row.c + 3 * (row.i + mx * (row.j + my * row.k));
    col_index = col.c + 3 * (col.i + mx * (col.j + my * col.k));
    PetscCall(DMDAGetAO(user->fda, &ao));
    PetscCall(AOApplicationToPetsc(ao, 1, &row_index));
    PetscCall(AOApplicationToPetsc(ao, 1, &col_index));
    PetscCall(DMCreateGlobalVector(user->fda, &column_basis));
    PetscCall(DMCreateGlobalVector(user->fda, &row_basis));
    PetscCall(DMCreateGlobalVector(user->fda, &product));
    PetscCall(VecSet(column_basis, 0.0)); PetscCall(VecSet(row_basis, 0.0));
    PetscCall(VecSetValue(column_basis, col_index, 1.0, INSERT_VALUES));
    PetscCall(VecSetValue(row_basis, row_index, 1.0, INSERT_VALUES));
    PetscCall(VecAssemblyBegin(column_basis)); PetscCall(VecAssemblyEnd(column_basis));
    PetscCall(VecAssemblyBegin(row_basis)); PetscCall(VecAssemblyEnd(row_basis));
    PetscCall(MatMult(preconditioning_matrix, column_basis, product));
    PetscCall(VecDot(row_basis, product, value));
    PetscCall(VecDestroy(&product)); PetscCall(VecDestroy(&row_basis));
    PetscCall(VecDestroy(&column_basis));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies the exact AIJ layout and preallocation derived from row classes. */
static PetscErrorCode AssertExactPointBlockMatrixAllocation(
    UserCtx *user, Mat matrix, PetscBool require_offrank_periodic)
{
    DMDALocalInfo info;
    MatInfo matrix_info;
    PetscMPIInt comm_size = 1;
    PetscInt matrix_rows, matrix_cols, local_rows, local_cols;
    PetscInt vector_size, vector_local_size, block_size;
    PetscInt ownership_start, ownership_end;
    PetscInt expected_local = 0, expected_global = 0;
    PetscInt offrank_periodic_local = 0, offrank_periodic_global = 0;
    PetscBool is_seq_aij = PETSC_FALSE, is_mpi_aij = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(DMDAGetLocalInfo(user->fda, &info));
    PetscCall(MatGetSize(matrix, &matrix_rows, &matrix_cols));
    PetscCall(MatGetLocalSize(matrix, &local_rows, &local_cols));
    PetscCall(VecGetSize(user->Ucont, &vector_size));
    PetscCall(VecGetLocalSize(user->Ucont, &vector_local_size));
    PetscCall(VecGetOwnershipRange(user->Ucont, &ownership_start, &ownership_end));
    PetscCall(MatGetBlockSize(matrix, &block_size));
    PetscCall(PicurvAssertIntEqual(vector_size, matrix_rows,
                                   "point-block global row dimension"));
    PetscCall(PicurvAssertIntEqual(vector_size, matrix_cols,
                                   "point-block global column dimension"));
    PetscCall(PicurvAssertIntEqual(vector_local_size, local_rows,
                                   "point-block local row dimension"));
    PetscCall(PicurvAssertIntEqual(vector_local_size, local_cols,
                                   "point-block local column dimension"));
    PetscCall(PicurvAssertIntEqual(3, block_size, "point-block logical block size"));
    PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)matrix), &comm_size));
    PetscCall(PetscObjectTypeCompare((PetscObject)matrix, MATSEQAIJ, &is_seq_aij));
    PetscCall(PetscObjectTypeCompare((PetscObject)matrix, MATMPIAIJ, &is_mpi_aij));
    PetscCall(PicurvAssertBool(
        comm_size == 1 ? is_seq_aij : is_mpi_aij,
        "point-block matrix must use the expected AIJ implementation"));
 
    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) {
                for (PetscInt component = 0; component < 3; ++component) {
                    PetscInt ri, rj, rk;
                    MomentumRowType type = ClassifyMomentumRow(
                        user, i, j, k, component, &ri, &rj, &rk);
                    expected_local += type == MOM_ROW_PHYSICAL ? 3 :
                                      type == MOM_ROW_PERIODIC_DUPLICATE ? 2 : 1;
                    if (type == MOM_ROW_PERIODIC_DUPLICATE) {
                        PetscInt representative;
                        PetscCall(TestStencilToGlobal(user,
                            (MatStencil){.i = ri, .j = rj, .k = rk, .c = component},
                            &representative));
                        if (representative < ownership_start || representative >= ownership_end)
                            ++offrank_periodic_local;
                    }
                }
            }
        }
    }
    PetscCallMPI(MPI_Allreduce(&expected_local, &expected_global, 1, MPIU_INT,
                               MPI_SUM, PetscObjectComm((PetscObject)matrix)));
    PetscCallMPI(MPI_Allreduce(&offrank_periodic_local, &offrank_periodic_global,
                               1, MPIU_INT, MPI_SUM,
                               PetscObjectComm((PetscObject)matrix)));
    PetscCall(MatGetInfo(matrix, MAT_GLOBAL_SUM, &matrix_info));
    PetscCall(PicurvAssertRealNear((PetscReal)expected_global,
        (PetscReal)matrix_info.nz_allocated, 0.0,
        "point-block matrix must allocate exactly the classified scalar pattern"));
    PetscCall(PicurvAssertRealNear((PetscReal)expected_global,
        (PetscReal)matrix_info.nz_used, 0.0,
        "point-block assembly must insert every classified structural entry"));
    PetscCall(PicurvAssertRealNear(0.0, (PetscReal)matrix_info.mallocs, 0.0,
        "point-block insertion must not reallocate matrix storage"));
    if (require_offrank_periodic && comm_size > 1)
        PetscCall(PicurvAssertBool((PetscBool)(offrank_periodic_global > 0),
            "MPI periodic fixture must exercise off-rank preallocation"));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Verifies the point-block model and common preconditioning-engine wiring.
 * @return PETSc error code.
 */
static PetscErrorCode TestPointBlockPreconditionerEngine(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    MomentumNewtonKrylovContext ctx = {0};
    MomentumPreconditionerDescription description = {
        MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN,
        MOM_NK_PC_STRUCTURE_POINT_BLOCK,
        0,
        0
    };
    MomentumPreconditionerEngine engine = {0};
    Mat preconditioning_matrix = NULL;
    Vec x = NULL, f = NULL;
    PetscInt block_size = 0, velocity_dof = 0;
    PetscReal matrix_norm = 0.0, reassembled_norm = 0.0, difference_norm = 0.0;
    MatStencil conditioned = {.i = 0, .j = 2, .k = 3, .c = 0};
    MatStencil homogeneous = {.i = 0, .j = 2, .k = 3, .c = 1};
    Mat saved_matrix = NULL, mffd = NULL;
    SNES snes = NULL;
    Vec direction = NULL, product = NULL, px = NULL;
    PetscReal px_norm = 0.0;
    KSP ksp = NULL;
    PC pc = NULL;
    Vec pc_rhs = NULL, pc_solution = NULL;
    MatInfo initial_allocation_info, repeated_allocation_info;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &f));
    PetscCall(VecCopy(user->Ucont, x));
    PetscCall(VecStrideSet(x, 0, 0.2));
    PetscCall(VecStrideSet(x, 1, 0.3));
    PetscCall(VecStrideSet(x, 2, 0.4));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f, &ctx));
    PetscCall(MomentumPreconditionerEngine_Create(user, NULL, &description, &engine));
    preconditioning_matrix = engine.preconditioning_matrix;
    PetscCall(PicurvAssertBool((PetscBool)(engine.model_ops == &frozen_momentum_point_block_ops),
                               "engine must select the frozen-momentum model callbacks"));
    PetscCall(PicurvAssertBool(engine.owns_preconditioning_matrix,
                               "point-block engine must own its separate matrix"));
    PetscCall(PicurvAssertBool((PetscBool)!engine.aliases_jacobian_operator,
                               "point-block matrix must not alias the Jacobian operator"));
    PetscCall(MomentumPreconditionerEngine_Assemble(&engine, user, x));
    PetscCall(AssertExactPointBlockMatrixAllocation(
        user, preconditioning_matrix, PETSC_FALSE));
    PetscCall(MatGetInfo(preconditioning_matrix, MAT_GLOBAL_SUM,
                         &initial_allocation_info));
    PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &matrix_norm));
    PetscCall(PicurvAssertBool((PetscBool)(matrix_norm > 0.0),
                               "model callback and common rows must insert matrix entries"));
    PetscCall(MatGetBlockSize(preconditioning_matrix, &block_size));
    PetscCall(PicurvAssertIntEqual(3, block_size, "point-block matrix block size"));
    PetscCall(DMDAGetInfo(user->fda, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                           &velocity_dof, NULL, NULL, NULL, NULL, NULL));
    PetscCall(PicurvAssertIntEqual(3, velocity_dof, "Newton velocity DMDA dof"));
 
    /* The coefficient formulas define columns by the differentiated Ucont
       component. Check an asymmetric pair so a row/column transpose cannot pass. */
    {
        const MatStencil row_i = {.i = 2, .j = 2, .k = 2, .c = 0};
        const MatStencil row_j = {.i = 2, .j = 2, .k = 2, .c = 1};
        MatStencil col_i = row_i, col_j = row_j;
        PetscReal ***aj = NULL;
        PetscReal AJip, AJjp;
        PetscScalar dFi_dUj = 0.0, dFj_dUi = 0.0;
 
        PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, &aj));
        AJip = 0.5 * (aj[2][2][2] + aj[2][2][3]);
        AJjp = 0.5 * (aj[2][2][2] + aj[2][3][2]);
        PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, &aj));
        PetscCall(PreconditionerMatrixStencilEntry(user, preconditioning_matrix,
            row_i, col_j, &dFi_dUj));
        PetscCall(PreconditionerMatrixStencilEntry(user, preconditioning_matrix,
            row_j, col_i, &dFj_dUi));
        PetscCall(PicurvAssertRealNear(0.5 * AJjp * 0.2,
            PetscRealPart(dFi_dUj), 1e-12,
            "row i, column j must differentiate with respect to U-j"));
        PetscCall(PicurvAssertRealNear(0.5 * AJip * 0.3,
            PetscRealPart(dFj_dUi), 1e-12,
            "row j, column i must differentiate with respect to U-i"));
    }
 
    /* Both fixed categories are exact identity rows, with no same-cell coupling. */
    for (PetscInt cc = 0; cc < 3; ++cc) {
        MatStencil col = conditioned; col.c = cc;
        PetscScalar conditioned_value = 0.0, homogeneous_value = 0.0;
        PetscCall(PreconditionerMatrixStencilEntry(user, preconditioning_matrix,
            conditioned, col, &conditioned_value));
        col = homogeneous; col.c = cc;
        PetscCall(PreconditionerMatrixStencilEntry(user, preconditioning_matrix,
            homogeneous, col, &homogeneous_value));
        PetscCall(PicurvAssertRealNear(cc == conditioned.c ? 1.0 : 0.0,
            PetscRealPart(conditioned_value), 1e-14, "conditioned row must be exact identity"));
        PetscCall(PicurvAssertRealNear(cc == homogeneous.c ? 1.0 : 0.0,
            PetscRealPart(homogeneous_value), 1e-14, "homogeneous row must be exact identity"));
    }
 
    /* A real residual and MFFD product cannot change subsequent assembly. */
    PetscCall(MatDuplicate(preconditioning_matrix, MAT_COPY_VALUES, &saved_matrix));
    PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    PetscCall(SNESSetDM(snes, user->fda));
    PetscCall(SNESSetFunction(snes, f, MomentumNewtonKrylov_FormResidual, &ctx));
    PetscCall(MatCreateSNESMF(snes, &mffd));
    PetscCall(VecDuplicate(x, &direction)); PetscCall(VecDuplicate(x, &product));
    PetscCall(VecSet(direction, 0.375));
    PetscCall(MomentumNewtonKrylov_FormResidual(snes, x, f, &ctx));
    PetscCall(MatMFFDSetBase(mffd, x, f));
    PetscCall(MatMult(mffd, direction, product));
    PetscCall(MomentumPreconditionerEngine_Assemble(&engine, user, x));
    PetscCall(MatAXPY(saved_matrix, -1.0, preconditioning_matrix, SAME_NONZERO_PATTERN));
    PetscCall(MatNorm(saved_matrix, NORM_FROBENIUS, &difference_norm));
    PetscCall(PicurvAssertRealNear(0.0, difference_norm, 1e-12,
        "assembly must be unchanged after residual and MFFD products"));
    PetscCall(VecDestroy(&product)); PetscCall(VecDestroy(&direction));
    PetscCall(MatDestroy(&mffd)); PetscCall(SNESDestroy(&snes)); PetscCall(MatDestroy(&saved_matrix));
 
    PetscCall(VecDuplicate(x, &px));
    PetscCall(MatMult(preconditioning_matrix, x, px));
    PetscCall(VecNorm(px, NORM_2, &px_norm));
    PetscCall(PicurvAssertBool((PetscBool)(!PetscIsInfOrNanReal(px_norm) && px_norm > 0.0),
        "point-block product must be finite and nonzero"));
    PetscCall(VecDestroy(&px));
 
    {
        PetscBool assembled = PETSC_FALSE;
        PetscCall(MatAssembled(preconditioning_matrix, &assembled));
        PetscCall(PicurvAssertBool(assembled, "engine must perform final matrix assembly"));
    }
    PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &matrix_norm));
    PetscCall(MatShift(preconditioning_matrix, 7.0));
    PetscCall(MomentumPreconditionerEngine_Assemble(&engine, user, x));
    PetscCall(MatGetInfo(preconditioning_matrix, MAT_GLOBAL_SUM,
                         &repeated_allocation_info));
    PetscCall(PicurvAssertRealNear((PetscReal)initial_allocation_info.nz_allocated,
        (PetscReal)repeated_allocation_info.nz_allocated, 0.0,
        "repeated assembly must retain exact allocated storage"));
    PetscCall(PicurvAssertRealNear((PetscReal)initial_allocation_info.mallocs,
        (PetscReal)repeated_allocation_info.mallocs, 0.0,
        "repeated assembly must not add insertion reallocations"));
    PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &reassembled_norm));
    PetscCall(PicurvAssertRealNear(matrix_norm, reassembled_norm, 1e-12,
                                   "repeated engine assembly must clear old entries"));
    {
        PetscScalar value = 0.0;
        PetscCall(PreconditionerMatrixStencilEntry(user, preconditioning_matrix,
            conditioned, conditioned, &value));
        PetscCall(PicurvAssertRealNear(1.0, PetscRealPart(value), 1e-14,
            "repeated engine assembly must clear old entries"));
    }
    PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    PetscCall(KSPSetOperators(ksp, preconditioning_matrix, preconditioning_matrix));
    PetscCall(KSPGetPC(ksp, &pc));
    PetscCall(MomentumPreconditionerEngine_ConfigurePetscPC(&engine, pc));
    PetscCall(MomentumPreconditionerEngine_ValidatePetscPC(&engine, pc));
    PetscCall(KSPSetUp(ksp));
    PetscCall(DMCreateGlobalVector(user->fda, &pc_rhs));
    PetscCall(DMCreateGlobalVector(user->fda, &pc_solution));
    PetscCall(VecSet(pc_rhs, 1.0));
    PetscCall(PCApply(pc, pc_rhs, pc_solution));
    PetscCall(VecDestroy(&pc_solution)); PetscCall(VecDestroy(&pc_rhs));
    PetscCall(KSPDestroy(&ksp));
    PetscCall(MomentumPreconditionerEngine_Destroy(&engine));
    PetscCall(PicurvAssertBool((PetscBool)(engine.preconditioning_matrix == NULL),
                               "engine destroy must clear its owned matrix"));
    PetscCall(VecDestroy(&f)); PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Proves that one periodic matrix row has only its exact +1/-1 pair. */
static PetscErrorCode AssertPeriodicPreconditionerRow(UserCtx *user, Mat matrix,
    MatStencil row, const char *message)
{
    AO ao = NULL;
    PetscInt ri, rj, rk, global_row, global_rep, lo, hi, local_checked = 0, checked = 0;
    PetscInt ncols = 0;
    const PetscInt *cols = NULL;
    const PetscScalar *values = NULL;
    MomentumRowType type;
 
    PetscFunctionBeginUser;
    type = ClassifyMomentumRow(user, row.i, row.j, row.k, row.c, &ri, &rj, &rk);
    PetscCall(PicurvAssertIntEqual(MOM_ROW_PERIODIC_DUPLICATE, type, message));
    global_row = row.c + 3 * (row.i + user->info.mx * (row.j + user->info.my * row.k));
    PetscCall(DMDAGetAO(user->fda, &ao));
    PetscCall(AOApplicationToPetsc(ao, 1, &global_row));
    PetscCall(MatGetOwnershipRange(matrix, &lo, &hi));
    if (global_row >= lo && global_row < hi) {
        PetscBool found_self = PETSC_FALSE, found_rep = PETSC_FALSE;
        PetscInt nonzero_entries = 0;
        PetscCall(TestStencilToGlobal(user,
            (MatStencil){.i = ri, .j = rj, .k = rk, .c = row.c}, &global_rep));
        PetscCall(MatGetRow(matrix, global_row, &ncols, &cols, &values));
        for (PetscInt n = 0; n < ncols; ++n) {
            if (PetscAbsScalar(values[n]) <= 1e-14) continue;
            nonzero_entries++;
            if (cols[n] == global_row) {
                found_self = PETSC_TRUE;
                PetscCall(PicurvAssertRealNear(1.0, PetscRealPart(values[n]), 1e-14,
                    "periodic row self entry"));
            } else if (cols[n] == global_rep) {
                found_rep = PETSC_TRUE;
                PetscCall(PicurvAssertRealNear(-1.0, PetscRealPart(values[n]), 1e-14,
                    "periodic row representative entry"));
            } else {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
                        "Periodic row contains an unintended nonzero column.");
            }
        }
        PetscCall(PicurvAssertIntEqual(2, nonzero_entries,
            "periodic row must contain exactly two numerical entries"));
        PetscCall(PicurvAssertBool((PetscBool)(found_self && found_rep), message));
        PetscCall(MatRestoreRow(matrix, global_row, &ncols, &cols, &values));
        local_checked = 1;
    }
    PetscCallMPI(MPI_Allreduce(&local_checked, &checked, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
    PetscCall(PicurvAssertIntEqual(1, checked, "periodic row must have exactly one matrix owner"));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Verifies Jacobian creation/registration and baseline alias ownership. */
static PetscErrorCode TestJacobianInterfaceAndBaselineAlias(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    MomentumNewtonKrylovContext ctx = {0};
    MomentumPreconditionerDescription description = {
        MOM_NK_PC_MODEL_NONE, MOM_NK_PC_STRUCTURE_NONE, 0, 0
    };
    SNES snes = NULL;
    Mat jacobian_operator = NULL, preconditioning_matrix = NULL;
    KSP ksp = NULL;
    PC pc = NULL;
    PetscInt rows = 0, cols = 0;
    const char *jacobian_prefix = NULL;
    const char *pc_type = NULL;
    PetscBool prefix_matches = PETSC_FALSE;
    PetscBool pc_is_none = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(fixed_wall_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    PetscCall(SNESSetDM(snes, user->fda));
    ctx.user = user;
    ctx.jacobian.type = MOM_NK_JACOBIAN_FINITE_DIFFERENCE;
    ctx.jacobian.finite_difference_mode = MOM_NK_FD_MODE_MATRIX_FREE;
    PetscCall(SNESSetFunction(snes, NULL, MomentumNewtonKrylov_FormResidual, &ctx));
    PetscCall(MomentumNewtonJacobian_Create(snes, &ctx.jacobian));
    PetscCall(MatGetOptionsPrefix(ctx.jacobian.jacobian_operator, &jacobian_prefix));
    PetscCall(PetscStrcmp(jacobian_prefix, "mom_nk_", &prefix_matches));
    PetscCall(PicurvAssertBool(prefix_matches,
                               "Jacobian interface must apply the application prefix"));
    PetscCall(MomentumPreconditionerEngine_Create(
        user, ctx.jacobian.jacobian_operator, &description, &ctx.preconditioning_engine));
    PetscCall(MomentumNewtonJacobian_Register(
        snes, &ctx.jacobian, &ctx.preconditioning_engine, &ctx));
    PetscCall(SNESGetKSP(snes, &ksp));
    PetscCall(KSPGetPC(ksp, &pc));
    PetscCall(MomentumPreconditionerEngine_ConfigurePetscPC(
        &ctx.preconditioning_engine, pc));
    PetscCall(MomentumPreconditionerEngine_ValidatePetscPC(
        &ctx.preconditioning_engine, pc));
    PetscCall(PCGetType(pc, &pc_type));
    PetscCall(PetscStrcmp(pc_type, PCNONE, &pc_is_none));
    PetscCall(PicurvAssertBool(pc_is_none,
                               "baseline engine must derive PETSc PCNONE"));
    PetscCall(SNESGetJacobian(snes, &jacobian_operator, &preconditioning_matrix, NULL, NULL));
    PetscCall(PicurvAssertBool(
        (PetscBool)(jacobian_operator == ctx.jacobian.jacobian_operator),
        "SNES must receive the Jacobian-interface MFFD operator"));
    PetscCall(PicurvAssertBool((PetscBool)(preconditioning_matrix == jacobian_operator),
                               "baseline preconditioning matrix must alias the Jacobian"));
    PetscCall(PicurvAssertBool(ctx.preconditioning_engine.aliases_jacobian_operator,
                               "baseline engine must record matrix aliasing"));
    PetscCall(PicurvAssertBool((PetscBool)!ctx.preconditioning_engine.owns_preconditioning_matrix,
                               "baseline engine must not own the Jacobian alias"));
    PetscCall(MomentumPreconditionerEngine_Destroy(&ctx.preconditioning_engine));
    PetscCall(MatGetSize(ctx.jacobian.jacobian_operator, &rows, &cols));
    PetscCall(PicurvAssertBool((PetscBool)(rows > 0 && cols > 0),
                               "destroying an alias engine must preserve the Jacobian"));
    PetscCall(MomentumNewtonJacobian_Destroy(&ctx.jacobian));
    PetscCall(SNESDestroy(&snes));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/** @brief Exercises engine-owned periodic duplicate rows on every MPI layout. */
static PetscErrorCode TestPointBlockPeriodicAssembly(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    MomentumNewtonKrylovContext ctx = {0};
    MomentumPreconditionerDescription description = {
        MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN,
        MOM_NK_PC_STRUCTURE_POINT_BLOCK,
        0,
        0
    };
    Vec x = NULL, f = NULL;
    PetscReal matrix_norm = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(BuildNewtonFixture(periodic_xyz_bcs, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &x));
    PetscCall(VecDuplicate(user->Ucont, &f));
    PetscCall(VecCopy(user->Ucont, x));
    ctx.user = user;
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, x, f, &ctx));
    PetscCall(MomentumPreconditionerEngine_Create(
        user, NULL, &description, &ctx.preconditioning_engine));
    PetscCall(MomentumPreconditionerEngine_Assemble(&ctx.preconditioning_engine, user, x));
    PetscCall(AssertExactPointBlockMatrixAllocation(
        user, ctx.preconditioning_engine.preconditioning_matrix, PETSC_TRUE));
    PetscCall(MatNorm(ctx.preconditioning_engine.preconditioning_matrix,
                      NORM_FROBENIUS, &matrix_norm));
    PetscCall(PicurvAssertBool((PetscBool)(matrix_norm > 0.0),
                               "periodic point-block assembly must produce a nonzero matrix"));
    PetscCall(AssertPeriodicPreconditionerRow(user,
        ctx.preconditioning_engine.preconditioning_matrix,
        (MatStencil){.i = 0, .j = 2, .k = 3, .c = 0},
        "single-axis periodic row must contain exact +1/-1 entries"));
    PetscCall(AssertPeriodicPreconditionerRow(user,
        ctx.preconditioning_engine.preconditioning_matrix,
        (MatStencil){.i = 0, .j = 0, .k = 3, .c = 1},
        "periodic intersection must contain exact +1/-1 entries"));
    PetscCall(AssertPeriodicPreconditionerRow(user,
        ctx.preconditioning_engine.preconditioning_matrix,
        (MatStencil){.i = 0, .j = 0, .k = 0, .c = 2},
        "periodic origin intersection must contain exact +1/-1 entries"));
    PetscCall(MomentumPreconditionerEngine_Destroy(&ctx.preconditioning_engine));
    PetscCall(VecDestroy(&f));
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(DestroyNewtonFixture(&simCtx, tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Runs the focused Newton--Krylov unit suite.
 * @param argc Command-line argument count.
 * @param argv Command-line argument vector.
 * @return Process exit status.
 */
int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"residual-repeatability-and-input-integrity", TestResidualRepeatabilityAndInputIntegrity},
        {"constraint-rows", TestConstraintRows},
        {"fixed-constraint-derivatives-all-faces", TestFixedConstraintDerivativesAllFaces},
        {"inlet-outlet-constraint-derivatives", TestInletOutletConstraintDerivatives},
        {"periodic-constraint-derivatives-and-intersections", TestPeriodicConstraintDerivativesAndIntersections},
        {"matrix-free-derivative", TestMatrixFreeDerivative},
        {"whole-operator-direct-jacobian", TestWholeOperatorDirectJacobian},
        {"periodic-operator-has-no-zero-rows", TestPeriodicOperatorHasNoZeroRows},
        {"zero-iteration-structured-logging", TestZeroIterationStructuredLogging},
        {"flat-channel-bdf1-startup", TestFlatChannelStartup},
        {"restart-and-continuation-solve", TestRestartAndContinuationSolve},
        {"small-solve-and-rollback", TestSmallSolveAndRollback},
        {"free-slip-walls-solve", TestFreeSlipWallsSolve},
        {"two-d-solve-holds-masked-component", TestTwoDSolveHoldsMaskedComponent},
        {"unsupported-configuration-fails-before-allocation", TestUnsupportedConfigurationFailsBeforeAllocation},
        {"post-allocation-failure-cleanup", TestPostAllocationFailureCleanup},
        {"linearization-config-parsing", TestLinearizationConfigParsing},
        {"point-block-preconditioner-engine", TestPointBlockPreconditionerEngine},
        {"jacobian-interface-and-baseline-alias", TestJacobianInterfaceAndBaselineAlias},
        {"point-block-periodic-assembly", TestPointBlockPeriodicAssembly},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv Newton Krylov tests");
    if (ierr) return (int)ierr;
    ierr = PicurvRunTests("unit-newton-krylov", cases, sizeof(cases) / sizeof(cases[0]));
    if (PetscFinalize()) return 1;
    return (int)ierr;
}