test_setup_lifecycle.c

Go to the documentation of this file.

/**
 * @file test_setup_lifecycle.c
 * @brief C unit tests for setup, initialization, and cleanup lifecycle entry points.
 */
 
#include "test_support.h"
 
#include "ParticleSwarm.h"
#include "ParticlePhysics.h"
#include "initialcondition.h"
#include "runloop.h"
#include "setup.h"
 
#include <stdio.h>
#include <string.h>
 
/**
 * @brief Asserts that a directory path exists and is readable.
 */
static PetscErrorCode AssertDirectoryExists(const char *path, const char *context)
{
    PetscBool exists = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(PetscTestDirectory(path, 'r', &exists));
    PetscCall(PicurvAssertBool(exists, context));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Writes one small temporary text file used by the partial-lifecycle context-only fixture.
 */
static PetscErrorCode WriteContextOnlyFile(const char *path, const char *contents)
{
    FILE *file = NULL;
 
    PetscFunctionBeginUser;
    file = fopen(path, "w");
    PetscCheck(file != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Failed to open '%s' for writing.", path);
    fputs(contents, file);
    fclose(file);
    PetscFunctionReturn(0);
}
 
/**
 * @brief Creates the control-file bundle needed for the context-only cleanup test.
 */
static PetscErrorCode PrepareContextOnlyConfig(char *tmpdir,
                                               size_t tmpdir_len,
                                               char *control_path,
                                               size_t control_path_len)
{
    char bcs_path[PETSC_MAX_PATH_LEN];
    char post_path[PETSC_MAX_PATH_LEN];
    char output_dir[PETSC_MAX_PATH_LEN];
    char log_dir[PETSC_MAX_PATH_LEN];
    char control_buffer[8192];
 
    PetscFunctionBeginUser;
    PetscCall(PicurvMakeTempDir(tmpdir, tmpdir_len));
    PetscCall(PetscSNPrintf(bcs_path, sizeof(bcs_path), "%s/bcs.run", tmpdir));
    PetscCall(PetscSNPrintf(post_path, sizeof(post_path), "%s/post.run", tmpdir));
    PetscCall(PetscSNPrintf(output_dir, sizeof(output_dir), "%s/results", tmpdir));
    PetscCall(PetscSNPrintf(log_dir, sizeof(log_dir), "%s/logs", tmpdir));
    PetscCall(PetscSNPrintf(control_path, control_path_len, "%s/test.control", tmpdir));
 
    PetscCall(WriteContextOnlyFile(
        bcs_path,
        "-Xi WALL noslip\n"
        "+Xi WALL noslip\n"
        "-Eta WALL noslip\n"
        "+Eta WALL noslip\n"
        "-Zeta INLET constant_velocity vx=0.0 vy=0.0 vz=1.5\n"
        "+Zeta OUTLET conservation\n"));
    PetscCall(WriteContextOnlyFile(
        post_path,
        "startTime = 0\n"
        "endTime = 1\n"
        "timeStep = 1\n"
        "output_particles = false\n"));
    PetscCall(PetscSNPrintf(
        control_buffer,
        sizeof(control_buffer),
        "-start_step 0\n"
        "-totalsteps 2\n"
        "-ren 100.0\n"
        "-dt 0.001\n"
        "-finit 1\n"
        "-ucont_x 0.0\n"
        "-ucont_y 0.0\n"
        "-ucont_z 1.5\n"
        "-bcs_files %s\n"
        "-profiling_timestep_mode off\n"
        "-profiling_final_summary true\n"
        "-postprocessing_config_file %s\n"
        "-grid\n"
        "-im 6\n"
        "-jm 6\n"
        "-km 6\n"
        "-xMins 0.0\n"
        "-xMaxs 1.0\n"
        "-yMins 0.0\n"
        "-yMaxs 1.0\n"
        "-zMins 0.0\n"
        "-zMaxs 1.0\n"
        "-rxs 1.0\n"
        "-rys 1.0\n"
        "-rzs 1.0\n"
        "-cgrids 0\n"
        "-nblk 1\n"
        "-euler_field_source solve\n"
        "-mom_solver_type DUALTIME_PICARD_RK4\n"
        "-mom_dt_rk4_residual_norm_noise_allowance_factor 1.07\n"
        "-mom_resid_atol 1.0e-8\n"
        "-mom_resid_rtol 1.0e-3\n"
        "-mom_nk_pic_monitor\n"
        "-mg_level 1\n"
        "-poisson 0\n"
        "-tio 0\n"
        "-particle_console_output_freq 0\n"
        "-logfreq 1\n"
        "-output_dir %s\n"
        "-restart_dir %s\n"
        /* This fixture owns its temporary tree and wants it wiped. That is an external
         * absolute path, so it authorizes the delete explicitly - the runtime guard
         * makes no exception for tests. */
        "-allow_unsafe_log_dir true\n"
        "-log_dir %s\n"
        "-numParticles 0\n"
        "-pinit 2\n",
        bcs_path,
        post_path,
        output_dir,
        output_dir,
        log_dir));
    PetscCall(WriteContextOnlyFile(control_path, control_buffer));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Finalizes and frees one lifecycle test context, then clears any PETSc options used to build it.
 */
static PetscErrorCode FreeLifecycleContext(SimCtx **simCtx_ptr)
{
    PetscFunctionBeginUser;
    PetscCall(PicurvDestroyRuntimeContext(simCtx_ptr));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Builds a full setup fixture through environment, grid, BC, and domain-rank initialization.
 */
static PetscErrorCode BuildLifecycleContext(PetscBool enable_particles, SimCtx **simCtx_out, char *tmpdir, size_t tmpdir_len)
{
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContext(NULL, enable_particles, simCtx_out, NULL, tmpdir, tmpdir_len));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Builds only the top-level simulation context used by partial-initialization cleanup tests.
 */
static PetscErrorCode BuildContextOnly(SimCtx **simCtx_out, char *tmpdir, size_t tmpdir_len)
{
    char control_path[PETSC_MAX_PATH_LEN];
    SimCtx *simCtx = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscOptionsClear(NULL));
    PetscCall(PrepareContextOnlyConfig(tmpdir, tmpdir_len, control_path, sizeof(control_path)));
    PetscCall(PetscOptionsSetValue(NULL, "-control_file", control_path));
    PetscCall(CreateSimulationContext(0, NULL, &simCtx));
    simCtx->exec_mode = EXEC_MODE_SOLVER;
    *simCtx_out = simCtx;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that setup refuses every flag whose feature is planned but not implemented.
 *
 * Multi-block coupling, immersed boundaries, moving bodies, moving frames and the
 * immersed-body flux corrections were all accepted before, and each ran a different
 * problem without saying so. The conductor refuses the YAML switches; this covers the
 * same flags arriving through a PETSc passthrough. Each flag is appended to an otherwise
 * valid control file, so a later duplicate overrides the fixture's own value.
 */
static PetscErrorCode TestSetupRejectsUnimplementedFeatureFlags(void)
{
    const char *refused[] = {"-nblk 2", "-imm 1", "-fsi 1", "-rfsi 1",
                             "-mframe 1", "-rframe 1", "-mhv 1", "-lv 1"};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(refused) / sizeof(refused[0]); ++n) {
        char tmpdir[PETSC_MAX_PATH_LEN];
        char control_path[PETSC_MAX_PATH_LEN];
        char message[128];
        SimCtx *simCtx = NULL;
        PetscErrorCode setup_ierr;
        FILE *file = NULL;
 
        PetscCall(PetscOptionsClear(NULL));
        PetscCall(PrepareContextOnlyConfig(tmpdir, sizeof(tmpdir), control_path, sizeof(control_path)));
        file = fopen(control_path, "a");
        PetscCheck(file != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Failed to reopen '%s'.", control_path);
        fprintf(file, "%s\n", refused[n]);
        fclose(file);
        PetscCall(PetscOptionsSetValue(NULL, "-control_file", control_path));
 
        PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
        setup_ierr = CreateSimulationContext(0, NULL, &simCtx);
        PetscCall(PetscPopErrorHandler());
        PetscCall(PetscSNPrintf(message, sizeof(message), "setup must refuse '%s'", refused[n]));
        PetscCall(PicurvAssertIntEqual(PETSC_ERR_SUP, setup_ierr, message));
 
        /* The refusal fires mid-parse, before any solver state exists; only the top-level
           context has been allocated. */
        PetscCall(PetscFree(simCtx));
        PetscCall(PetscOptionsClear(NULL));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Builds a context from the fixture control file with extra lines appended.
 * @param[in]  extra_lines Control-file lines appended after the fixture's own.
 * @param[out] simCtx_out  Context on success; the partially built context on failure.
 * @param[out] ierr_out    CreateSimulationContext's return code, errors suppressed.
 * @param[out] tmpdir      Fixture directory, removed by the caller.
 * @param[in]  tmpdir_len  Capacity of `tmpdir`.
 */
static PetscErrorCode TryContextWithExtraLines(const char *extra_lines, SimCtx **simCtx_out,
                                               PetscErrorCode *ierr_out, char *tmpdir,
                                               size_t tmpdir_len)
{
    char control_path[PETSC_MAX_PATH_LEN];
    FILE *file = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PetscOptionsClear(NULL));
    PetscCall(PrepareContextOnlyConfig(tmpdir, tmpdir_len, control_path, sizeof(control_path)));
    file = fopen(control_path, "a");
    PetscCheck(file != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Failed to reopen '%s'.", control_path);
    fputs(extra_lines, file);
    fclose(file);
    PetscCall(PetscOptionsSetValue(NULL, "-control_file", control_path));
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    *ierr_out = CreateSimulationContext(0, NULL, simCtx_out);
    PetscCall(PetscPopErrorHandler());
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the IEM mixing constant defaults to 2.0, follows -iem_constant,
 *        accepts zero (no micromixing), and refuses a negative or non-finite value.
 */
static PetscErrorCode TestIEMConstantIsConfigurable(void)
{
    const char *cases[] = {"", "-iem_constant 3.5\n", "-iem_constant 0\n", "-iem_constant -1\n",
                           "-iem_constant nan\n", "-iem_constant inf\n"};
    const PetscReal expected[] = {2.0, 3.5, 0.0, -1.0, -1.0, -1.0};  /* -1: refused */
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(cases) / sizeof(cases[0]); ++n) {
        char tmpdir[PETSC_MAX_PATH_LEN];
        SimCtx *simCtx = NULL;
        PetscErrorCode setup_ierr = 0;
 
        PetscCall(TryContextWithExtraLines(cases[n], &simCtx, &setup_ierr, tmpdir, sizeof(tmpdir)));
        if (expected[n] >= 0.0) {
            PetscCall(PicurvAssertIntEqual(0, setup_ierr, "a non-negative -iem_constant must be accepted"));
            PetscCall(PicurvAssertRealNear(expected[n], simCtx->iem_constant, 0.0,
                                           "the IEM constant must follow -iem_constant, default 2.0"));
        } else {
            PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_OUTOFRANGE, setup_ierr,
                                           "a negative or non-finite -iem_constant must be refused"));
        }
        PetscCall(FreeLifecycleContext(&simCtx));
        PetscCall(PetscOptionsClear(NULL));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Ensures a run that solves no pressure builds one grid level whatever -mg_level asks,
 *        and a solved flow keeps the levels it asks for.
 */
static PetscErrorCode TestUnsolvedFlowUsesOneGridLevel(void)
{
    const char *cases[] = {"-mg_level 3\n", "-mg_level 3\n-euler_field_source analytical\n",
                           "-mg_level 3\n-euler_field_source load\n"};
    const PetscInt expected[] = {3, 1, 1};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(cases) / sizeof(cases[0]); ++n) {
        char tmpdir[PETSC_MAX_PATH_LEN];
        SimCtx *simCtx = NULL;
        PetscErrorCode setup_ierr = 0;
 
        PetscCall(TryContextWithExtraLines(cases[n], &simCtx, &setup_ierr, tmpdir, sizeof(tmpdir)));
        PetscCall(PicurvAssertIntEqual(0, setup_ierr, "the multigrid level options must be accepted"));
        PetscCall(PicurvAssertIntEqual(expected[n], simCtx->mglevels,
                                       "only a solved flow keeps its coarse grid levels"));
        PetscCall(FreeLifecycleContext(&simCtx));
        PetscCall(PetscOptionsClear(NULL));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Ensures setup refuses a particle count whose widest swarm field PetscInt cannot index.
 * @details Particle fields have up to three components, so the limit is PETSC_MAX_INT / 3
 *          particles; with 64-bit indices no representable -numParticles reaches it.
 */
static PetscErrorCode TestParticleCountMustFitPetscInt(void)
{
    char tmpdir[PETSC_MAX_PATH_LEN], line[64];
    SimCtx *simCtx = NULL;
    PetscErrorCode setup_ierr = 0;
 
    PetscFunctionBeginUser;
    if (sizeof(PetscInt) > 4) PetscFunctionReturn(0);
    PetscCall(PetscSNPrintf(line, sizeof(line), "-numParticles %" PetscInt_FMT "\n", (PetscInt)(PETSC_MAX_INT / 3) + 1));
    PetscCall(TryContextWithExtraLines(line, &simCtx, &setup_ierr, tmpdir, sizeof(tmpdir)));
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_OUTOFRANGE, setup_ierr,
                                   "a particle count past PETSC_MAX_INT / 3 must be refused with 32-bit indices"));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscCall(PetscOptionsClear(NULL));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Ensures setup propagates invalid post intervals before either executable can loop.
 */
static PetscErrorCode TestSetupRejectsNonpositivePostInterval(void)
{
    const char *options[] = {"-timeStep 0\n", "-timeStep -1\n"};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(options) / sizeof(options[0]); ++n) {
        SimCtx *simCtx = NULL;
        char tmpdir[PETSC_MAX_PATH_LEN];
        PetscErrorCode setup_ierr;
        PetscCall(TryContextWithExtraLines(options[n], &simCtx, &setup_ierr, tmpdir, sizeof(tmpdir)));
        PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_OUTOFRANGE, setup_ierr,
                    "setup must propagate rejection of a nonpositive post interval"));
        PetscCall(FreeLifecycleContext(&simCtx));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Checks a configured particle initial value through the production setup path.
 * @details `Psi = x + 2*pid` is read from the options, compiled, and applied once positions
 *          are final; it varies across particles wherever they are seeded. Particles fill the
 *          volume, so cells next to every face are occupied. Every particle must carry
 *          exactly its value, and the t=0 scatter must already reflect it: the
 *          count-weighted Eulerian `Psi` sums to the particles' total, so the first IEM
 *          update relaxes toward the right mean, and each non-periodic face's dummy cells
 *          carry the adjacent cell's mean.
 */
static PetscErrorCode TestConfiguredParticleInitialValue(void)
{
    const char *options = "-numParticles 400\n"
                          "-pinit 1\n"
                          "-particle_fields_count 1\n"
                          "-particle_fields_0_name Psi\n"
                          "-particle_fields_0_expr_0 \"x + 2*pid\"\n";
    SimCtx          *simCtx = NULL;
    UserCtx         *user = NULL;
    char             tmpdir[PETSC_MAX_PATH_LEN];
    PetscInt         nlocal = 0;
    const PetscReal  *psi = NULL, *position = NULL;
    const PetscInt64 *pid = NULL;
    PetscReal         particle_sum = 0.0, eulerian_sum = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(NULL, PETSC_TRUE, options, &simCtx, &user,
                                                       tmpdir, sizeof(tmpdir)));
    PetscCall(PicurvAssertBool((PetscBool)(simCtx->particleFieldPlan != NULL), "the plan is read from the options"));
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(InitializeParticleSwarm(simCtx));
    PetscCall(PerformInitializedParticleSetup(simCtx));
 
    PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
    PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&position));
    PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
    for (PetscInt p = 0; p < nlocal; ++p) {
        const PetscReal expected = position[3 * p] * simCtx->scaling.L_ref + 2.0 * (PetscReal)pid[p];
        PetscCall(PicurvAssertRealNear(expected, psi[p], 1.0e-12, "each particle carries its initial value"));
        particle_sum += psi[p];
    }
    PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
    PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&position));
    PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &particle_sum, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
    PetscCall(VecDot(user->Psi, user->ParticleCount, &eulerian_sum));
    PetscCall(PicurvAssertRealNear(particle_sum, eulerian_sum, 1.0e-9,
                                   "the t=0 scatter carries the initial value to the Eulerian mean"));
    {
        /* The scatter leaves empty dummy cells at zero; the production path must then give
           each non-periodic face's dummy the adjacent cell's mean. */
        const DMDALocalInfo info = user->info;
        const PetscReal   ***mean = NULL, ***count = NULL;
        PetscReal           checked = 0.0;
 
        PetscCall(DMDAVecGetArrayRead(user->da, user->Psi, &mean));
        PetscCall(DMDAVecGetArrayRead(user->da, user->ParticleCount, &count));
        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++) {
            /* A flat-face dummy cell and its adjacent interior cell, one axis at a time. */
            const PetscInt  index[3] = {i, j, k}, extent[3] = {info.mx, info.my, info.mz};
            PetscInt        dummy_axes = 0, axis = -1;
 
            for (PetscInt d = 0; d < 3; d++) {
                if (index[d] == 0 || index[d] == extent[d] - 1) { dummy_axes++; axis = d; }
            }
            if (dummy_axes != 1) continue;
            if (user->boundary_faces[2 * axis].mathematical_type == PERIODIC) continue;
            {
                const PetscInt step = (index[axis] == 0) ? 1 : -1;
                const PetscInt ai = i + (axis == 0 ? step : 0), aj = j + (axis == 1 ? step : 0),
                               ak = k + (axis == 2 ? step : 0);
 
                if (count[ak][aj][ai] < 0.5) continue;
                PetscCall(PicurvAssertRealNear(mean[ak][aj][ai], mean[k][j][i], 0.0,
                                               "the scatter fills each dummy cell with the adjacent cell mean"));
                checked += 1.0;
            }
        }
        PetscCall(DMDAVecRestoreArrayRead(user->da, user->ParticleCount, &count));
        PetscCall(DMDAVecRestoreArrayRead(user->da, user->Psi, &mean));
        PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, &checked, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
        PetscCall(PicurvAssertBool((PetscBool)(checked > 0.0),
                                   "the fixture must seed particles next to a non-periodic face"));
    }
 
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Checks configured IEM relaxation through the production swarm update.
 * @details Each particle must relax toward the mean of the particles now in its cell, at the
 *          configured rate. A zero constant switches micromixing off, and must leave Psi
 *          bit-identical: alternating values 0.1 and 0.7 do not survive `m + (0.1 - m) * 1`
 *          exactly in mixed cells, so the exact comparison shows the update is skipped, not
 *          run at a zero rate.
 */
static PetscErrorCode TestConfiguredIEMUpdatesSwarm(void)
{
    const char *options[] = {"", "-iem_constant 3.5\n", "-iem_constant 0\n"};
    const PetscReal constants[] = {2.0, 3.5, 0.0};
 
    PetscFunctionBeginUser;
    for (size_t n = 0; n < sizeof(constants) / sizeof(constants[0]); ++n) {
        SimCtx *simCtx = NULL;
        UserCtx *user = NULL;
        char tmpdir[PETSC_MAX_PATH_LEN];
        PetscInt nlocal;
        PetscReal *psi = NULL, *diffusivity = NULL, *psi_initial = NULL, *cell_mean = NULL;
        PetscInt *cell = NULL;
        const char *psi_name = ParticleFieldName(PARTICLE_FIELD_ID_PSI);
        const char *diff_name = ParticleFieldName(PARTICLE_FIELD_ID_DIFFUSIVITY);
        const char *cell_name = ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID);
 
        PetscCall(PicurvBuildTinyRuntimeContextWithOptions(NULL, PETSC_TRUE, options[n],
                    &simCtx, &user, tmpdir, sizeof(tmpdir)));
        PetscCall(InitializeEulerianState(simCtx));
        PetscCall(InitializeParticleSwarm(simCtx));
        PetscCall(PerformInitializedParticleSetup(simCtx));
        simCtx->dt = 0.5;
        /* Controlled nonzero test data; no public scalar initializer is implied. A unit
         * Jacobian makes the mixing length 1, and alternating particle values give every
         * cell with two or more particles a mean the particles must relax toward. */
        PetscCall(VecSet(user->lAj, 1.0));
        PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
        PetscCall(PicurvAssertBool((PetscBool)(nlocal > 0), "IEM fixture must contain particles"));
        PetscCall(DMSwarmGetField(user->swarm, psi_name, NULL, NULL, (void **)&psi));
        PetscCall(DMSwarmGetField(user->swarm, diff_name, NULL, NULL, (void **)&diffusivity));
        PetscCall(DMSwarmGetField(user->swarm, cell_name, NULL, NULL, (void **)&cell));
        PetscCall(PetscMalloc2(nlocal, &psi_initial, nlocal, &cell_mean));
        for (PetscInt p = 0; p < nlocal; ++p) {
            psi[p] = psi_initial[p] = (p % 2) ? 0.7 : 0.1;
            diffusivity[p] = 0.2;
        }
        /* Each particle's own-cell mean, from the particles where they are now. */
        for (PetscInt p = 0; p < nlocal; ++p) {
            PetscReal sum = 0.0, count = 0.0;
            for (PetscInt q = 0; q < nlocal; ++q) {
                if (cell[3*q] == cell[3*p] && cell[3*q+1] == cell[3*p+1] && cell[3*q+2] == cell[3*p+2]) {
                    sum += psi_initial[q]; count += 1.0;
                }
            }
            cell_mean[p] = sum / count;
        }
        PetscCall(DMSwarmRestoreField(user->swarm, cell_name, NULL, NULL, (void **)&cell));
        PetscCall(DMSwarmRestoreField(user->swarm, diff_name, NULL, NULL, (void **)&diffusivity));
        PetscCall(DMSwarmRestoreField(user->swarm, psi_name, NULL, NULL, (void **)&psi));
        PetscCall(UpdateAllParticleFields(user));
        PetscCall(DMSwarmGetField(user->swarm, psi_name, NULL, NULL, (void **)&psi));
        for (PetscInt p = 0; p < nlocal; ++p) {
            if (constants[n] == 0.0) {
                PetscCall(PicurvAssertBool((PetscBool)(psi[p] == psi_initial[p]),
                            "a zero IEM constant must leave Psi bit-identical"));
            } else {
                PetscCall(PicurvAssertRealNear(cell_mean[p] + (psi_initial[p] - cell_mean[p]) * PetscExpReal(-constants[n] * 0.2 * 0.5),
                            psi[p], 1.e-12, "configured IEM constant must relax each particle toward its cell's current mean"));
            }
        }
        PetscCall(DMSwarmRestoreField(user->swarm, psi_name, NULL, NULL, (void **)&psi));
        PetscCall(PetscFree2(psi_initial, cell_mean));
        PetscCall(FreeLifecycleContext(&simCtx));
        PetscCall(PicurvRemoveTempDir(tmpdir));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Checks that IEM relaxes each particle toward its own cell's scattered mean.
 * @details Particles start at random 0 or 1, so cell means differ from cell to cell and
 *          from the zero in the boundary storage slots. Relaxing toward the particle's own
 *          cell mean at a rate uniform within the cell conserves every cell's scalar sum,
 *          so the total is unchanged while the spread shrinks; reading a neighbour's mean,
 *          or a boundary slot, breaks the conservation.
 */
static PetscErrorCode TestIEMRelaxesTowardOwnCellMean(void)
{
    const char *options = "-particle_fields_count 1\n"
                          "-particle_fields_0_name Psi\n"
                          "-particle_fields_0_expr_0 \"where(uniform() < 0.5, 0, 1)\"\n";
    SimCtx     *simCtx = NULL;
    UserCtx    *user = NULL;
    char        tmpdir[PETSC_MAX_PATH_LEN];
    PetscInt    nlocal = 0;
    PetscReal  *psi = NULL, *diffusivity = NULL;
    PetscReal   sums[2][2] = {{0.0, 0.0}, {0.0, 0.0}}; /* [before|after][sum|sum of squares] */
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(NULL, PETSC_TRUE, options, &simCtx, &user,
                                                       tmpdir, sizeof(tmpdir)));
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(InitializeParticleSwarm(simCtx));
    PetscCall(PerformInitializedParticleSetup(simCtx));
    simCtx->dt = 0.5;
 
    PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
    PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    PetscCall(DMSwarmGetField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
    for (PetscInt p = 0; p < nlocal; ++p) {
        diffusivity[p] = 0.2;
        sums[0][0] += psi[p];
        sums[0][1] += psi[p] * psi[p];
    }
    PetscCall(DMSwarmRestoreField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
    PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
 
    PetscCall(UpdateAllParticleFields(user));
 
    PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    for (PetscInt p = 0; p < nlocal; ++p) {
        sums[1][0] += psi[p];
        sums[1][1] += psi[p] * psi[p];
    }
    PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
    PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, sums, 4, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 
    PetscCall(PicurvAssertBool((PetscBool)(sums[0][0] > 0.0 && sums[0][0] < (PetscReal)simCtx->np),
                               "the fixture must start with both values present"));
    PetscCall(PicurvAssertRealNear(sums[0][0], sums[1][0], 1.0e-9,
                                   "relaxing toward each particle's own cell mean conserves the scalar"));
    PetscCall(PicurvAssertBool((PetscBool)(sums[1][1] < sums[0][1] - 1.0e-6),
                               "the update must reduce the scalar's spread"));
 
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Checks that post-processing loads the particle averages a checkpoint holds.
 * @details A case with `restart_mode: init` reseeds particles on a solver restart, so the
 *          solver must not restore the averages of the particles it discards; the
 *          postprocessor reads the same control settings but visualizes what was saved, so it
 *          must load them, and refresh the ghosted copy that nodal averaging reads.
 */
static PetscErrorCode TestPostprocessorLoadsSavedParticleAverages(void)
{
    const char *options = "-particle_fields_count 1\n"
                          "-particle_fields_0_name Psi\n"
                          "-particle_fields_0_expr_0 \"where(uniform() < 0.5, 0, 1)\"\n"
                          "-particle_restart_mode init\n";
    SimCtx     *simCtx = NULL;
    UserCtx    *user = NULL;
    char        tmpdir[PETSC_MAX_PATH_LEN];
    Vec         saved = NULL, local_saved = NULL;
    PetscReal   saved_norm = 0.0, difference = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContextWithOptions(NULL, PETSC_TRUE, options, &simCtx, &user,
                                                       tmpdir, sizeof(tmpdir)));
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(InitializeParticleSwarm(simCtx));
    PetscCall(PerformInitializedParticleSetup(simCtx));
    PetscCall(VecDuplicate(user->Psi, &saved));
    PetscCall(VecCopy(user->Psi, saved));
    PetscCall(VecNorm(saved, NORM_INFINITY, &saved_norm));
    PetscCall(PicurvAssertBool((PetscBool)(saved_norm > 0.0), "the fixture must scatter a nonzero Psi"));
    PetscCall(WriteCheckpointBundle(simCtx, "test"));
 
    PetscCall(PetscCalloc1(1, &simCtx->pps));
    PetscCall(PetscStrncpy(simCtx->pps->source_dir, simCtx->output_dir, sizeof(simCtx->pps->source_dir)));
    simCtx->exec_mode = EXEC_MODE_POSTPROCESSOR;
    PetscCall(VecSet(user->Psi, 0.0));
    PetscCall(VecSet(user->lPsi, 0.0));
    PetscCall(ReadSimulationFields(user, simCtx->step));
 
    PetscCall(VecAXPY(saved, -1.0, user->Psi));
    PetscCall(VecNorm(saved, NORM_INFINITY, &difference));
    PetscCall(PicurvAssertRealNear(0.0, difference, 1.0e-12, "post-processing restores the saved Psi average"));
    PetscCall(DMCreateLocalVector(user->da, &local_saved));
    PetscCall(DMGlobalToLocalBegin(user->da, user->Psi, INSERT_VALUES, local_saved));
    PetscCall(DMGlobalToLocalEnd(user->da, user->Psi, INSERT_VALUES, local_saved));
    PetscCall(VecAXPY(local_saved, -1.0, user->lPsi));
    PetscCall(VecNorm(local_saved, NORM_INFINITY, &difference));
    PetscCall(PicurvAssertRealNear(0.0, difference, 1.0e-12, "the ghosted Psi is refreshed for nodal averaging"));
 
    simCtx->exec_mode = EXEC_MODE_SOLVER;
    PetscCall(PetscFree(simCtx->pps));
    PetscCall(VecDestroy(&local_saved));
    PetscCall(VecDestroy(&saved));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that the shared richer runtime fixture mirrors normalized production setup contracts.
 */
static PetscErrorCode TestSharedRuntimeFixtureContracts(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PetscCall(PicurvBuildTinyRuntimeContext(NULL, PETSC_FALSE, &simCtx, &user, tmpdir, sizeof(tmpdir)));
 
    PetscCall(PicurvAssertBool((PetscBool)(user->bbox.min_coords.x >= -2.0e-6), "runtime fixture bbox xmin should stay inside normalized domain tolerance"));
    PetscCall(PicurvAssertBool((PetscBool)(user->bbox.max_coords.x <= 1.0 + 2.0e-6), "runtime fixture bbox xmax should stay inside normalized domain tolerance"));
    PetscCall(PicurvAssertBool((PetscBool)(user->bbox.max_coords.y <= 1.0 + 2.0e-6), "runtime fixture bbox ymax should stay inside normalized domain tolerance"));
    PetscCall(PicurvAssertBool((PetscBool)(user->bbox.max_coords.z <= 1.0 + 2.0e-6), "runtime fixture bbox zmax should stay inside normalized domain tolerance"));
    PetscCall(PicurvAssertBool((PetscBool)(simCtx->bboxlist != NULL), "runtime fixture should gather bboxlist through SetupDomainRankInfo"));
    PetscCall(PicurvAssertBool((PetscBool)(user->RankCellInfoMap != NULL), "runtime fixture should gather rank-cell ownership metadata"));
    PetscCall(PicurvAssertBool((PetscBool)(simCtx->bboxlist[simCtx->rank].max_coords.x >= user->bbox.max_coords.x - 1.0e-10),
                               "runtime fixture rank bbox entry should include the local bbox extent"));
 
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Asserts that one block's stored local info matches what its DMs report.
 * @details Grid, rank and particle code read user->info in place of DMDAGetLocalInfo on
 *          user->da and user->fda, which relies on both DMs sharing one decomposition.
 */
static PetscErrorCode AssertStoredLocalInfoMatchesDMs(UserCtx *user)
{
    DMDALocalInfo dm_info[2];
    const char   *names[2] = {"da", "fda"};
 
    PetscFunctionBeginUser;
    PetscCall(DMDAGetLocalInfo(user->da, &dm_info[0]));
    PetscCall(DMDAGetLocalInfo(user->fda, &dm_info[1]));
    for (PetscInt d = 0; d < 2; ++d) {
        const DMDALocalInfo *a = &dm_info[d], *b = &user->info;
        const PetscInt got[15]  = {a->mx, a->my, a->mz, a->xs, a->ys, a->zs, a->xm, a->ym, a->zm,
                                   a->gxs, a->gys, a->gzs, a->gxm, a->gym, a->gzm};
        const PetscInt want[15] = {b->mx, b->my, b->mz, b->xs, b->ys, b->zs, b->xm, b->ym, b->zm,
                                   b->gxs, b->gys, b->gzs, b->gxm, b->gym, b->gzm};
        char msg[128];
        for (PetscInt f = 0; f < 15; ++f) {
            PetscCall(PetscSNPrintf(msg, sizeof(msg), "user->info range %" PetscInt_FMT " should match %s", f, names[d]));
            PetscCall(PicurvAssertIntEqual(want[f], got[f], msg));
        }
        PetscCall(PicurvAssertBool((PetscBool)(a->bx == b->bx && a->by == b->by && a->bz == b->bz),
                                   "user->info boundary types should match the DM"));
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests the core setup lifecycle through environment, grid, BC, rank-info, and Eulerian-state initialization.
 */
static PetscErrorCode TestSetupLifecycleCoreSolverSetup(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
    char results_dir[PETSC_MAX_PATH_LEN];
    char logs_dir[PETSC_MAX_PATH_LEN];
    Cmpnts ***ucont = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(BuildLifecycleContext(PETSC_FALSE, &simCtx, tmpdir, sizeof(tmpdir)));
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    PetscCall(PetscSNPrintf(results_dir, sizeof(results_dir), "%s/results", tmpdir));
    PetscCall(PetscSNPrintf(logs_dir, sizeof(logs_dir), "%s/logs", tmpdir));
    PetscCall(PicurvAssertBool((PetscBool)(user->da != NULL), "SetupGridAndSolvers should allocate da"));
    PetscCall(PicurvAssertBool((PetscBool)(user->fda != NULL), "SetupGridAndSolvers should allocate coordinate DM"));
    PetscCall(PicurvAssertIntEqual(user->IM + 1, user->info.mx, "DM node count should match IM+1"));
    for (PetscInt level = 0; level < simCtx->usermg.mglevels; ++level) {
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            PetscCall(AssertStoredLocalInfoMatchesDMs(&simCtx->usermg.mgctx[level].user[bi]));
        }
    }
    PetscCall(PicurvAssertBool(user->inletFaceDefined, "SetupBoundaryConditions should identify the inlet face"));
    PetscCall(PicurvAssertIntEqual(BC_FACE_NEG_Z, user->identifiedInletBCFace, "NEG_Z should be the configured inlet face"));
    PetscCall(PicurvAssertBool((PetscBool)(simCtx->bboxlist != NULL), "SetupDomainRankInfo should build bboxlist"));
    PetscCall(PicurvAssertBool((PetscBool)(user->RankCellInfoMap != NULL), "SetupDomainRankInfo should build rank-cell decomposition map"));
    PetscCall(AssertDirectoryExists(results_dir, "SetupSimulationEnvironment should create the output directory"));
    PetscCall(AssertDirectoryExists(logs_dir, "SetupSimulationEnvironment should create the log directory"));
 
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
    PetscCall(PicurvAssertBool((PetscBool)(ucont[1][1][1].z > 0.0),
                               "InitializeEulerianState should seed a positive inlet-aligned interior field"));
    PetscCall(PicurvAssertRealNear(ucont[1][1][1].z, ucont[2][2][2].z, 1.0e-12,
                                   "InitializeEulerianState should initialize a spatially consistent interior field"));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
 
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests particle-swarm initialization and deterministic settlement on a tiny fully initialized case.
 */
static PetscErrorCode TestSetupLifecycleParticleInitialization(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
    PetscInt nlocal = 0;
    PetscInt *cell_ids = NULL;
    PetscInt *status = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(BuildLifecycleContext(PETSC_TRUE, &simCtx, tmpdir, sizeof(tmpdir)));
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(InitializeParticleSwarm(simCtx));
    PetscCall(PerformInitializedParticleSetup(simCtx));
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    PetscCall(PicurvAssertBool((PetscBool)(user->swarm != NULL), "InitializeParticleSwarm should allocate the solver swarm"));
    PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
    PetscCall(PicurvAssertIntEqual(8, nlocal, "single-rank lifecycle particle setup should own all seeded particles"));
 
    PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
    PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
    for (PetscInt p = 0; p < nlocal; ++p) {
        PetscCall(PicurvAssertBool((PetscBool)(cell_ids[3 * p + 0] >= 0), "settled particles should have a valid i cell id"));
        PetscCall(PicurvAssertBool((PetscBool)(cell_ids[3 * p + 1] >= 0), "settled particles should have a valid j cell id"));
        PetscCall(PicurvAssertBool((PetscBool)(cell_ids[3 * p + 2] >= 0), "settled particles should have a valid k cell id"));
        PetscCall(PicurvAssertIntEqual(ACTIVE_AND_LOCATED, status[p], "settled particles should be marked ACTIVE_AND_LOCATED"));
    }
    PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
    PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
 
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests that step-zero scalar verification emits scatter_metrics.csv during initialized setup.
 */
static PetscErrorCode TestSetupLifecycleScatterMetricsAtStepZero(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN];
    char metrics_path[PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PetscCall(BuildLifecycleContext(PETSC_TRUE, &simCtx, tmpdir, sizeof(tmpdir)));
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    PetscCall(PetscStrncpy(simCtx->eulerianSource, "analytical", sizeof(simCtx->eulerianSource)));
    PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
    simCtx->StepsToRun = 0;
    simCtx->StartStep = 0;
    simCtx->verificationScalar.enabled = PETSC_TRUE;
    PetscCall(PetscStrncpy(simCtx->verificationScalar.mode,
                           "analytical",
                           sizeof(simCtx->verificationScalar.mode)));
    PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
                           "CONSTANT",
                           sizeof(simCtx->verificationScalar.profile)));
    simCtx->verificationScalar.value = 1.25;
 
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(InitializeParticleSwarm(simCtx));
    PetscCall(PerformInitializedParticleSetup(simCtx));
 
    PetscCall(PetscSNPrintf(metrics_path, sizeof(metrics_path), "%s/results/analysis/metrics/scatter_metrics.csv", tmpdir));
    PetscCall(PicurvAssertFileExists(metrics_path,
                                     "PerformInitializedParticleSetup should emit scatter_metrics.csv for a zero-step scalar verification run"));
    PetscCall(PicurvAssertBool((PetscBool)(user->swarm != NULL),
                               "zero-step scalar verification setup should leave the swarm initialized"));
 
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscCall(FreeLifecycleContext(&simCtx));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests standalone RNG initialization helpers and minimal-context cleanup.
 */
static PetscErrorCode TestSetupLifecycleRandomGeneratorsAndCleanup(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    PetscRandom randx = NULL, randy = NULL, randz = NULL;
    PetscRandom rand_i = NULL, rand_j = NULL, rand_k = NULL;
    PetscScalar sample_x = 0.0, sample_i = 0.0;
    PetscReal seconds = 0.0;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
    user->bbox.min_coords.x = 0.0;
    user->bbox.min_coords.y = 0.0;
    user->bbox.min_coords.z = 0.0;
    user->bbox.max_coords.x = 1.0;
    user->bbox.max_coords.y = 2.0;
    user->bbox.max_coords.z = 3.0;
 
    PetscCall(InitializeRandomGenerators(user, &randx, &randy, &randz));
    PetscCall(InitializeLogicalSpaceRNGs(simCtx->particleRandomSeed, &rand_i, &rand_j, &rand_k));
    PetscCall(PicurvAssertBool(RuntimeWalltimeGuardParsePositiveSeconds("12.5", &seconds), "RuntimeWalltimeGuardParsePositiveSeconds should parse positive numeric strings"));
    PetscCall(PicurvAssertRealNear(12.5, seconds, 1.0e-12, "RuntimeWalltimeGuardParsePositiveSeconds parsed value"));
 
    PetscCall(PetscRandomGetValue(randx, &sample_x));
    PetscCall(PetscRandomGetValue(rand_i, &sample_i));
    PetscCall(PicurvAssertBool((PetscBool)(PetscRealPart(sample_x) >= 0.0 && PetscRealPart(sample_x) <= 1.0),
                               "InitializeRandomGenerators should honor the configured bbox interval"));
    PetscCall(PicurvAssertBool((PetscBool)(PetscRealPart(sample_i) >= 0.0 && PetscRealPart(sample_i) <= 1.0),
                               "InitializeLogicalSpaceRNGs should generate logical coordinates in [0,1]"));
 
    PetscCall(PetscRandomDestroy(&randx));
    PetscCall(PetscRandomDestroy(&randy));
    PetscCall(PetscRandomDestroy(&randz));
    PetscCall(PetscRandomDestroy(&rand_i));
    PetscCall(PetscRandomDestroy(&rand_j));
    PetscCall(PetscRandomDestroy(&rand_k));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Tests cleanup after partial and fuller setup states without requiring unsupported double-finalization behavior.
 */
static PetscErrorCode TestSetupLifecycleCleanupAcrossInitializationStates(void)
{
    SimCtx *context_only = NULL;
    SimCtx *grid_only = NULL;
    char context_tmpdir[PETSC_MAX_PATH_LEN];
    char grid_tmpdir[PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PetscCall(BuildContextOnly(&context_only, context_tmpdir, sizeof(context_tmpdir)));
    PetscCall(PicurvAssertBool((PetscBool)(context_only != NULL), "CreateSimulationContext should allocate the top-level SimCtx"));
    PetscCall(PicurvAssertIntEqual(1, context_only->block_number, "CreateSimulationContext should parse the configured block count"));
    PetscCall(PicurvAssertIntEqual(MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK, context_only->mom_solver_type,
                                   "deprecated RK4 selector should normalize to the Jameson solver enum"));
    PetscCall(PicurvAssertRealNear(1.07, context_only->mom_dt_jameson_residual_norm_noise_allowance_factor, 1.0e-12,
                                   "deprecated RK4 residual-noise option should populate the Jameson control"));
    PetscCall(PicurvAssertRealNear(1.0e-8, context_only->mom_resid_atol, 1.0e-12,
                                   "momentum absolute residual tolerance should be parsed"));
    PetscCall(PicurvAssertRealNear(1.0e-3, context_only->mom_resid_rtol, 1.0e-12,
                                   "momentum relative residual tolerance should be parsed"));
    PetscCall(PicurvAssertBool(context_only->mom_nk_monitor_history,
                               "Newton structured history monitor should be parsed"));
    PetscCall(PetscOptionsClear(NULL));
    PetscCall(PicurvRemoveTempDir(context_tmpdir));
    PetscCall(FreeLifecycleContext(&context_only));
 
    PetscCall(BuildContextOnly(&grid_only, grid_tmpdir, sizeof(grid_tmpdir)));
    PetscCall(SetupSimulationEnvironment(grid_only));
    PetscCall(SetupGridAndSolvers(grid_only));
    PetscCall(PicurvAssertBool((PetscBool)(grid_only->usermg.mgctx[grid_only->usermg.mglevels - 1].user->Ucont != NULL),
                               "SetupGridAndSolvers should allocate baseline Eulerian vectors"));
    PetscCall(PicurvRemoveTempDir(grid_tmpdir));
    PetscCall(FreeLifecycleContext(&grid_only));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Verifies catalog completeness, stable name resolution, and runtime view binding.
 */
static PetscErrorCode TestFieldCatalogMetadataAndViews(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    FieldId parsed_id = FIELD_ID_INVALID;
    FieldView view;
    const FieldDescriptor *descriptor = NULL;
    Vec coordinates = NULL;
    Vec local_coordinates = NULL;
    PetscErrorCode unknown_name_ierr;
    PetscErrorCode unavailable_view_ierr;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
 
    for (PetscInt raw_id = 0; raw_id < FIELD_ID_COUNT; ++raw_id) {
        PetscCall(FieldGetDescriptor((FieldId)raw_id, &descriptor));
        PetscCall(PicurvAssertIntEqual(raw_id, (PetscInt)descriptor->id,
                                       "catalog entry should retain its declared FieldId"));
        PetscCall(PicurvAssertBool((PetscBool)(descriptor->canonical_name != NULL &&
                                               descriptor->canonical_name[0] != '\0'),
                                   "every field descriptor should have a canonical name"));
        PetscCall(PicurvAssertBool((PetscBool)(descriptor->dof == 1 ||
                                               descriptor->dof == 2 ||
                                               descriptor->dof == 3),
                                   "every ghost-updatable catalog field should declare a supported dof"));
        PetscCall(PicurvAssertBool((PetscBool)((descriptor->capabilities &
                                               FIELD_CAPABILITY_GHOST_UPDATE) != 0u),
                                   "every phase-one field should advertise ghost-update capability"));
        PetscCall(FieldIdFromName(descriptor->canonical_name, &parsed_id));
        PetscCall(PicurvAssertIntEqual(raw_id, (PetscInt)parsed_id,
                                       "canonical field names should round-trip to their IDs"));
        PetscCall(PicurvAssertBool((PetscBool)(strcmp(FieldCanonicalName((FieldId)raw_id),
                                                     descriptor->canonical_name) == 0),
                                   "FieldCanonicalName should return catalog-owned canonical text"));
    }
 
    PetscCall(FieldIdFromName("eddy viscosity", &parsed_id));
    PetscCall(PicurvAssertIntEqual(FIELD_ID_NU_T, parsed_id,
                                   "Nu_t compatibility alias should resolve case-insensitively"));
    PetscCall(FieldIdFromName("cs", &parsed_id));
    PetscCall(PicurvAssertIntEqual(FIELD_ID_CS, parsed_id,
                                   "CS compatibility alias should resolve case-insensitively"));
    PetscCall(FieldIdFromName("X-Face-Centers", &parsed_id));
    PetscCall(PicurvAssertIntEqual(FIELD_ID_CENTX, parsed_id,
                                   "diagnostic face-center name should resolve to the canonical field"));
 
    PetscCall(FieldGetDescriptor(FIELD_ID_UCONT, &descriptor));
    PetscCall(PicurvAssertIntEqual(FIELD_LAYOUT_COMPONENT_STAGGERED, descriptor->layout,
                                   "Ucont should retain component-staggered topology metadata"));
    PetscCall(PicurvAssertIntEqual(FIELD_SYNC_COMPONENT_STAGGERED, descriptor->sync_class,
                                   "Ucont should retain component-staggered ghost repair metadata"));
    PetscCall(FieldGetDescriptor(FIELD_ID_CSI, &descriptor));
    PetscCall(PicurvAssertIntEqual(FIELD_LAYOUT_I_FACE, descriptor->layout,
                                   "Csi should retain I-face topology metadata"));
    PetscCall(PicurvAssertIntEqual(FIELD_SYNC_I_FACE, descriptor->sync_class,
                                   "Csi should retain I-face ghost repair metadata"));
    PetscCall(FieldGetDescriptor(FIELD_ID_IAJ, &descriptor));
    PetscCall(PicurvAssertIntEqual(1, descriptor->dof,
                                   "IAj should retain scalar face-field metadata"));
 
    PetscCall(FieldGetView(user, FIELD_ID_P, &view));
    PetscCall(PicurvAssertBool((PetscBool)(view.dm == user->da &&
                                           view.global_vec == user->P &&
                                           view.local_vec == user->lP),
                               "P field view should bind the existing scalar DM and Vec pair"));
    PetscCall(FieldGetView(user, FIELD_ID_UCAT, &view));
    PetscCall(PicurvAssertBool((PetscBool)(view.dm == user->fda &&
                                           view.global_vec == user->Ucat &&
                                           view.local_vec == user->lUcat),
                               "Ucat field view should bind the existing vector DM and Vec pair"));
 
    PetscCall(DMGetCoordinates(user->da, &coordinates));
    PetscCall(DMGetCoordinatesLocal(user->da, &local_coordinates));
    PetscCall(FieldGetView(user, FIELD_ID_COORDINATES, &view));
    PetscCall(PicurvAssertBool((PetscBool)(view.dm == user->fda &&
                                           view.global_vec == coordinates &&
                                           view.local_vec == local_coordinates),
                               "coordinate field view should bind PETSc-owned coordinate vectors"));
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    unknown_name_ierr = FieldIdFromName("NotARegisteredEulerianField", &parsed_id);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_UNKNOWN_TYPE, unknown_name_ierr,
                                   "unknown field names should fail at the name-ingress boundary"));
    PetscCall(PicurvAssertIntEqual(FIELD_ID_INVALID, parsed_id,
                                   "failed field-name resolution should return FIELD_ID_INVALID"));
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    unavailable_view_ierr = FieldGetView(user, FIELD_ID_CS, &view);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_WRONGSTATE, unavailable_view_ierr,
                                   "catalogued optional fields without storage should report unavailable state"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Verifies the separate DMSwarm field catalog and its Eulerian bridge metadata.
 */
static PetscErrorCode TestParticleFieldCatalogMetadata(void)
{
    ParticleFieldId parsed_id = PARTICLE_FIELD_ID_INVALID;
    const ParticleFieldDescriptor *descriptor = NULL;
    PetscErrorCode unknown_name_ierr;
 
    PetscFunctionBeginUser;
    for (PetscInt raw_id = 0; raw_id < PARTICLE_FIELD_ID_COUNT; ++raw_id) {
        PetscCall(ParticleFieldGetDescriptor((ParticleFieldId)raw_id, &descriptor));
        PetscCall(PicurvAssertIntEqual(raw_id, (PetscInt)descriptor->id,
                                       "particle catalog entry should retain its declared ID"));
        PetscCall(PicurvAssertBool((PetscBool)(descriptor->canonical_name != NULL &&
                                               descriptor->canonical_name[0] != '\0'),
                                   "every particle field should have a canonical DMSwarm name"));
        PetscCall(PicurvAssertBool((PetscBool)(descriptor->components == 1 ||
                                               descriptor->components == 3),
                                   "persistent particle fields should declare one or three components"));
        PetscCall(ParticleFieldIdFromName(descriptor->canonical_name, &parsed_id));
        PetscCall(PicurvAssertIntEqual(raw_id, (PetscInt)parsed_id,
                                       "canonical particle field names should round-trip to their IDs"));
        PetscCall(PicurvAssertBool((PetscBool)(strcmp(ParticleFieldName((ParticleFieldId)raw_id),
                                                     descriptor->canonical_name) == 0),
                                   "ParticleFieldName should return catalog-owned PETSc text"));
    }
 
    PetscCall(ParticleFieldIdFromName("particlevelocity", &parsed_id));
    PetscCall(PicurvAssertIntEqual(PARTICLE_FIELD_ID_VELOCITY, parsed_id,
                                   "particle velocity alias should resolve case-insensitively"));
    PetscCall(ParticleFieldIdFromName("Migration Status", &parsed_id));
    PetscCall(PicurvAssertIntEqual(PARTICLE_FIELD_ID_LOCATION_STATUS, parsed_id,
                                   "migration-status output alias should resolve to the PETSc field"));
 
    PetscCall(ParticleFieldGetDescriptor(PARTICLE_FIELD_ID_PSI, &descriptor));
    PetscCall(PicurvAssertBool((PetscBool)((descriptor->capabilities &
                                           PARTICLE_FIELD_CAPABILITY_EULERIAN_SCATTER) != 0u),
                               "Psi should advertise its active particle-to-Eulerian scatter"));
    PetscCall(PicurvAssertIntEqual(FIELD_ID_PSI, descriptor->eulerian_scatter_target,
                                   "particle Psi should target Eulerian Psi"));
    PetscCall(ParticleFieldGetDescriptor(PARTICLE_FIELD_ID_PID, &descriptor));
    PetscCall(PicurvAssertIntEqual(PARTICLE_FIELD_REGISTRATION_PETSC, descriptor->registration,
                                   "particle IDs should be marked as PETSc-managed"));
    PetscCall(PicurvAssertIntEqual(PETSC_INT64, descriptor->data_type,
                                   "particle IDs should retain 64-bit storage metadata"));
 
    PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
    unknown_name_ierr = ParticleFieldIdFromName("NotARegisteredParticleField", &parsed_id);
    PetscCall(PetscPopErrorHandler());
    PetscCall(PicurvAssertIntEqual(PETSC_ERR_ARG_UNKNOWN_TYPE, unknown_name_ierr,
                                   "unknown particle names should fail at the text-ingress boundary"));
    PetscCall(PicurvAssertIntEqual(PARTICLE_FIELD_ID_INVALID, parsed_id,
                                   "failed particle name resolution should return the invalid ID"));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Verifies typed ghost updates preserve scalar, vector, face, and staggered behavior.
 */
static PetscErrorCode TestFieldCatalogTypedGhostUpdates(void)
{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
 
    PetscCall(VecSet(user->P, 4.25));
    PetscCall(VecSet(user->lP, -1.0));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_P));
    PetscCall(PicurvAssertVecConstant(user->lP, 4.25, 1.0e-12,
                                      "typed scalar ghost update should scatter P"));
 
    PetscCall(VecSet(user->Ucat, 7.5));
    PetscCall(VecSet(user->lUcat, -1.0));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCAT));
    PetscCall(PicurvAssertVecConstant(user->lUcat, 7.5, 1.0e-12,
                                      "typed cell-vector ghost update should scatter Ucat"));
 
    PetscCall(VecSet(user->IAj, 2.0));
    PetscCall(VecSet(user->lIAj, -1.0));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_IAJ));
    PetscCall(PicurvAssertVecConstant(user->lIAj, 2.0, 1.0e-12,
                                      "typed scalar face-field ghost update should scatter IAj"));
 
    PetscCall(VecSet(user->Ucont, 3.0));
    PetscCall(VecSet(user->lUcont, -1.0));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCONT));
    PetscCall(PicurvAssertVecConstant(user->lUcont, 3.0, 1.0e-12,
                                      "typed component-staggered ghost update should scatter Ucont"));
 
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Runs the unit-setup PETSc test binary.
 */
int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"setup-lifecycle-core-solver-setup", TestSetupLifecycleCoreSolverSetup},
        {"setup-lifecycle-particle-initialization", TestSetupLifecycleParticleInitialization},
        {"setup-lifecycle-scatter-metrics-step-zero", TestSetupLifecycleScatterMetricsAtStepZero},
        {"setup-lifecycle-random-generators-and-cleanup", TestSetupLifecycleRandomGeneratorsAndCleanup},
        {"setup-lifecycle-cleanup-across-initialization-states", TestSetupLifecycleCleanupAcrossInitializationStates},
        {"setup-rejects-unimplemented-feature-flags", TestSetupRejectsUnimplementedFeatureFlags},
        {"iem-constant-is-configurable", TestIEMConstantIsConfigurable},
        {"unsolved-flow-uses-one-grid-level", TestUnsolvedFlowUsesOneGridLevel},
        {"particle-count-must-fit-petscint", TestParticleCountMustFitPetscInt},
        {"setup-rejects-nonpositive-post-interval", TestSetupRejectsNonpositivePostInterval},
        {"configured-iem-updates-swarm", TestConfiguredIEMUpdatesSwarm},
        {"configured-particle-initial-value", TestConfiguredParticleInitialValue},
        {"iem-relaxes-toward-own-cell-mean", TestIEMRelaxesTowardOwnCellMean},
        {"postprocessor-loads-saved-particle-averages", TestPostprocessorLoadsSavedParticleAverages},
        {"shared-runtime-fixture-contracts", TestSharedRuntimeFixtureContracts},
        {"field-catalog-metadata-and-views", TestFieldCatalogMetadataAndViews},
        {"particle-field-catalog-metadata", TestParticleFieldCatalogMetadata},
        {"field-catalog-typed-ghost-updates", TestFieldCatalogTypedGhostUpdates},
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv setup lifecycle tests");
    if (ierr) {
        return (int)ierr;
    }
 
    ierr = PicurvRunTests("unit-setup", cases, sizeof(cases) / sizeof(cases[0]));
    if (ierr) {
        PetscFinalize();
        return (int)ierr;
    }
 
    ierr = PetscFinalize();
    return (int)ierr;
}