setup.c

Go to the documentation of this file.

/**
 * @file setup.c  //  Setup code for running any simulation 
 * @brief Test program for DMSwarm interpolation using the fdf-curvIB method.
 * Provides the setup to start any simulation with DMSwarm and DMDAs.
 **/
 
#include <ctype.h>
#include <errno.h>
 
#include "setup.h"
#include "statistics_config.h"
#include "statistics_accumulator.h"
#include "ParticleInitialConditions.h"
#include "verification_sources.h"
#include <unistd.h>
#include <limits.h>
#include <stdlib.h>
 
#define PICURV_STRINGIZE_(x) #x
#define PICURV_STRINGIZE(x) PICURV_STRINGIZE_(x)
#if defined(PETSC_USE_DEBUG)
#define PICURV_PETSC_MODE "debug"
#else
#define PICURV_PETSC_MODE "optimized"
#endif
#ifndef PETSC_ARCH
#define PETSC_ARCH ""
#endif
#ifndef PETSC_DIR
#define PETSC_DIR ""
#endif
#define PICURV_PETSC_STAMP_MARKER "PICURV_PETSC_BUILD:"
 
/* The PETSc this binary was compiled against. It is printed by --version and also read
 * from the file by the conductor, which needs it most when a library mismatch keeps the
 * binary from starting at all. */
static const char picurv_petsc_build_stamp[] =
    PICURV_PETSC_STAMP_MARKER "petsc "
    PICURV_STRINGIZE(PETSC_VERSION_MAJOR) "." PICURV_STRINGIZE(PETSC_VERSION_MINOR) "."
    PICURV_STRINGIZE(PETSC_VERSION_SUBMINOR) " " PICURV_PETSC_MODE " " PETSC_ARCH " " PETSC_DIR;
 
/**
 * @brief Implementation of \ref PicurvHandleVersionArgument().
 * @details The public argument and return contract is documented in `include/setup.h`.
 * @see PicurvHandleVersionArgument()
 */
int PicurvHandleVersionArgument(int argc, char **argv, const char *executable_name)
{
    for (int index = 1; index < argc; ++index) {
        if (!strcmp(argv[index], "--version") || !strcmp(argv[index], "-version")) {
            printf("%s %s+g%.12s%s\n", executable_name, PICURV_RELEASE_VERSION,
                   PICURV_GIT_COMMIT,
                   !strcmp(PICURV_BUILD_DIRTY, "true") ? ".dirty" : "");
            printf("%s\n", picurv_petsc_build_stamp + sizeof(PICURV_PETSC_STAMP_MARKER) - 1);
            return 1;
        }
    }
    return 0;
}
 
/**
 * @brief Implementation of \ref RuntimeWalltimeGuardParsePositiveSeconds().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see RuntimeWalltimeGuardParsePositiveSeconds()
 */
PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
{
    char   *endptr = NULL;
    double parsed_value;
 
    if (seconds_out) *seconds_out = 0.0;
    if (!text || text[0] == '\0') return PETSC_FALSE;
 
    errno        = 0;
    parsed_value = strtod(text, &endptr);
    if (endptr == text || errno == ERANGE || !isfinite(parsed_value) || parsed_value <= 0.0) {
        return PETSC_FALSE;
    }
 
    while (*endptr != '\0' && isspace((unsigned char)*endptr)) {
        endptr++;
    }
    if (*endptr != '\0') return PETSC_FALSE;
 
    if (seconds_out) *seconds_out = (PetscReal)parsed_value;
    return PETSC_TRUE;
}
 
/**
 * @brief Implementation of \ref InitializeSolutionConvergenceState().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see InitializeSolutionConvergenceState()
 */
PetscErrorCode InitializeSolutionConvergenceState(SimCtx *simCtx)
{
    UserCtx *user = NULL;
    PetscInt history_capacity = 0;
 
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL.");
    if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
    if (!simCtx->solutionConvergenceEnabled) PetscFunctionReturn(0);
    if (!simCtx->usermg.mgctx) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "Multigrid hierarchy must exist before initializing solution convergence storage.");
    }
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    if (!user) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "Finest-level UserCtx must exist before initializing solution convergence storage.");
    }
 
    simCtx->solutionConvergenceSamplesRecorded = 0;
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC) {
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
                                   &user[bi].solutionConvergencePeriodicUcatRef));
            PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
                                   &user[bi].solutionConvergencePeriodicPRef));
            for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                PetscCall(VecDuplicate(user[bi].Ucat, &user[bi].solutionConvergencePeriodicUcatRef[phase]));
                PetscCall(VecSet(user[bi].solutionConvergencePeriodicUcatRef[phase], 0.0));
                PetscCall(VecDuplicate(user[bi].P, &user[bi].solutionConvergencePeriodicPRef[phase]));
                PetscCall(VecSet(user[bi].solutionConvergencePeriodicPRef[phase], 0.0));
            }
        }
    }
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY) {
        history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
        PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanSpeedHistory));
        PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanKEHistory));
    }
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref DestroySolutionConvergenceState().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see DestroySolutionConvergenceState()
 */
PetscErrorCode DestroySolutionConvergenceState(SimCtx *simCtx)
{
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx) PetscFunctionReturn(0);
    if (simCtx->fieldStatisticsWindows) {
        PetscCall(PetscFree(simCtx->fieldStatisticsWindows));
        simCtx->fieldStatisticsWindows = NULL;
        simCtx->fieldStatisticsWindowCount = 0;
    }
    if (!simCtx->usermg.mgctx) PetscFunctionReturn(0);
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    if (user) {
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            if (user[bi].solutionConvergencePeriodicUcatRef) {
                for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                    if (user[bi].solutionConvergencePeriodicUcatRef[phase]) {
                        PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicUcatRef[phase]));
                    }
                }
                PetscCall(PetscFree(user[bi].solutionConvergencePeriodicUcatRef));
                user[bi].solutionConvergencePeriodicUcatRef = NULL;
            }
            if (user[bi].solutionConvergencePeriodicPRef) {
                for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                    if (user[bi].solutionConvergencePeriodicPRef[phase]) {
                        PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicPRef[phase]));
                    }
                }
                PetscCall(PetscFree(user[bi].solutionConvergencePeriodicPRef));
                user[bi].solutionConvergencePeriodicPRef = NULL;
            }
        }
    }
 
    if (simCtx->solutionConvergenceMeanSpeedHistory) {
        PetscCall(PetscFree(simCtx->solutionConvergenceMeanSpeedHistory));
        simCtx->solutionConvergenceMeanSpeedHistory = NULL;
    }
    if (simCtx->solutionConvergenceMeanKEHistory) {
        PetscCall(PetscFree(simCtx->solutionConvergenceMeanKEHistory));
        simCtx->solutionConvergenceMeanKEHistory = NULL;
    }
    simCtx->solutionConvergenceSamplesRecorded = 0;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LESConfigSetDefaults"
/**
 * @brief Implementation of \ref LESConfigSetDefaults().
 * @details Full API contract is documented with the header declaration in
 *          `include/setup.h`.
 * @see LESConfigSetDefaults()
 */
PetscErrorCode LESConfigSetDefaults(LESConfig *config)
{
    PetscFunctionBeginUser;
    PetscCheck(config != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "LES configuration destination cannot be NULL.");
 
    config->dynamic_frequency       = 1;
    config->constant_cs             = 0.03;
    // Vreman (2004) recommends c = 2.5 Cs^2, 0.07 for Cs = 0.17; Nicoud & Ducros (1999)
    // calibrate C_w near 0.5-0.6 against the same isotropic decay. Both are calibrations
    // from isotropic turbulence, not universal values.
    config->vreman_coefficient      = 0.07;
    config->wale_coefficient        = 0.5;
    config->filter_width_model      = LES_FILTER_WIDTH_CUBE_ROOT_VOLUME;
    config->test_filter_kernel      = LES_TEST_FILTER_VOLUME_WEIGHTED_BOX;
    config->test_filter_width_ratio = 2.0;
    // Local averaging makes no homogeneity claim, so it is the only default that is
    // correct on every geometry. Homogeneous and global averaging are opted into.
    config->averaging_mode          = LES_AVERAGING_LOCAL;
    config->averaging_direction[0]  = PETSC_FALSE;
    config->averaging_direction[1]  = PETSC_FALSE;
    config->averaging_direction[2]  = PETSC_FALSE;
    config->clip_mode               = LES_CLIP_CLAMP;
    // Roughly twice the physical Smagorinsky constant: high enough to catch a
    // diverging coefficient, low enough that it does not shape the distribution.
    config->max_cs                  = 0.3;
    config->min_viscosity_ratio     = 0.0;
    config->yoshizawa_ci            = 0.09;
    config->diagnostics_enabled     = PETSC_FALSE;
    config->diagnostics_cadence     = 1;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ParseLESConfiguration"
/**
 * @brief Reads every LES closure parameter from the generated control file.
 *
 * Grouped into one routine because the parameters constrain one another: the
 * averaging directions only mean something under homogeneous averaging, and the
 * clip ceiling only under clamping. Selectors arrive as integers that the Python
 * layer has already normalized from their user-facing names, which is the same
 * contract `-les`, `-pinit`, and `-interpolation_method` follow.
 *
 * Values that would make the closure ill-defined are rejected here as well as in
 * Python, so a hand-written control file fails with a reason rather than producing
 * a plausible-looking coefficient.
 *
 * @param[in,out] simCtx Simulation context whose `les_config` is populated.
 * @return PetscErrorCode 0 on success, or `PETSC_ERR_ARG_OUTOFRANGE` for a value
 *         outside its admissible range.
 */
static PetscErrorCode ParseLESConfiguration(SimCtx *simCtx)
{
    LESConfig *config = &simCtx->les_config;
    PetscInt   selector;
    char       directions[8] = "";
    PetscBool  found = PETSC_FALSE;
 
    PetscFunctionBeginUser;
 
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_constant_cs", &config->constant_cs, NULL));
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_vreman_coefficient", &config->vreman_coefficient, NULL));
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_wale_coefficient", &config->wale_coefficient, NULL));
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_dynamic_frequency", &config->dynamic_frequency, NULL));
 
    selector = (PetscInt)config->filter_width_model;
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_filter_width", &selector, NULL));
    PetscCheck(selector >= LES_FILTER_WIDTH_CUBE_ROOT_VOLUME && selector <= LES_FILTER_WIDTH_SCOTTI,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_filter_width must be 0 (cube_root_volume), 1 (geometric_mean), 2 (max_edge), or 3 (scotti); received %" PetscInt_FMT ".",
               selector);
    config->filter_width_model = (LESFilterWidthModel)selector;
 
    selector = (PetscInt)config->test_filter_kernel;
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_test_filter_kernel", &selector, NULL));
    PetscCheck(selector >= LES_TEST_FILTER_VOLUME_WEIGHTED_BOX && selector <= LES_TEST_FILTER_SIMPSON_IK,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_test_filter_kernel must be 0 (volume_weighted_box) or 1 (simpson_ik); received %" PetscInt_FMT ".",
               selector);
    config->test_filter_kernel = (LESTestFilterKernel)selector;
 
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_test_filter_width_ratio",
                                  &config->test_filter_width_ratio, NULL));
    PetscCheck(config->test_filter_width_ratio > 1.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_test_filter_width_ratio must exceed 1; a test filter no wider than the grid "
               "filter leaves the dynamic procedure nothing to measure. Received %g.",
               (double)config->test_filter_width_ratio);
 
    selector = (PetscInt)config->averaging_mode;
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_averaging_mode", &selector, NULL));
    PetscCheck(selector >= LES_AVERAGING_LOCAL && selector <= LES_AVERAGING_GLOBAL,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_averaging_mode must be 0 (local), 1 (homogeneous), or 2 (global); received %" PetscInt_FMT ".",
               selector);
    config->averaging_mode = (LESAveragingMode)selector;
 
    // A subset of "ijk". Left empty under homogeneous averaging, the periodic axes
    // are used instead, which is what makes the common cases need no configuration.
    PetscCall(PetscOptionsGetString(NULL, NULL, "-les_averaging_directions",
                                    directions, sizeof(directions), &found));
    if (found) {
        for (size_t index = 0; directions[index] != '\0'; index++) {
            switch (directions[index]) {
            case 'i': config->averaging_direction[0] = PETSC_TRUE; break;
            case 'j': config->averaging_direction[1] = PETSC_TRUE; break;
            case 'k': config->averaging_direction[2] = PETSC_TRUE; break;
            default:
                SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                        "-les_averaging_directions accepts only the characters i, j, and k; "
                        "received '%s'.", directions);
            }
        }
    }
 
    selector = (PetscInt)config->clip_mode;
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_clip_mode", &selector, NULL));
    PetscCheck(selector >= LES_CLIP_CLAMP && selector <= LES_CLIP_NONE,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_clip_mode must be 0 (clamp), 1 (clip_negative), or 2 (none); received %" PetscInt_FMT ".",
               selector);
    config->clip_mode = (LESClipMode)selector;
 
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_clip_max_cs", &config->max_cs, NULL));
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_min_viscosity_ratio",
                                  &config->min_viscosity_ratio, NULL));
    PetscCall(PetscOptionsGetReal(NULL, NULL, "-les_yoshizawa_ci", &config->yoshizawa_ci, NULL));
    PetscCall(PetscOptionsGetBool(NULL, NULL, "-les_diagnostics",
                                  &config->diagnostics_enabled, NULL));
    PetscCall(PetscOptionsGetInt(NULL, NULL, "-les_diagnostics_cadence",
                                 &config->diagnostics_cadence, NULL));
 
    PetscCheck(config->dynamic_frequency > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_dynamic_frequency must be positive; received %" PetscInt_FMT ".",
               config->dynamic_frequency);
    PetscCheck(config->constant_cs >= 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_constant_cs must be nonnegative; received %g.", (double)config->constant_cs);
    PetscCheck(config->vreman_coefficient >= 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_vreman_coefficient must be nonnegative; received %g.",
               (double)config->vreman_coefficient);
    PetscCheck(config->wale_coefficient >= 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_wale_coefficient must be nonnegative; received %g.",
               (double)config->wale_coefficient);
    PetscCheck(config->max_cs >= 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_clip_max_cs must be nonnegative; received %g.", (double)config->max_cs);
    PetscCheck(config->min_viscosity_ratio >= 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_min_viscosity_ratio must be nonnegative; received %g.",
               (double)config->min_viscosity_ratio);
    PetscCheck(config->diagnostics_cadence > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-les_diagnostics_cadence must be positive; received %" PetscInt_FMT ".",
               config->diagnostics_cadence);
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "CreateSimulationContext"
 
/**
 * @brief Implementation of \ref CreateSimulationContext().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see CreateSimulationContext()
 */
PetscErrorCode CreateSimulationContext(int argc, char **argv, SimCtx **p_simCtx)
{
    PetscErrorCode ierr;
    (void)argc;
    (void)argv;
    SimCtx *simCtx;
    char control_filename[PETSC_MAX_PATH_LEN] = ""; // Temporary placeholder for control file name.
    PetscBool control_flg; // Temporary placeholder for control file tag existence check flag.
    PetscBool particle_console_output_freq_flg = PETSC_FALSE;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // === 1. Allocate the Context Struct and Set ALL Defaults ==================
    ierr = PetscNew(p_simCtx); CHKERRQ(ierr);
    simCtx = *p_simCtx;
 
    // --- Group 1: Parallelism & MPI ---
    simCtx->rank = 0; simCtx->size = 1;
 
    // --- Group 2: Simulation Control, Time, and I/O ---
    simCtx->step = 0; simCtx->ti = 0.0; simCtx->StartStep = 0; simCtx->StepsToRun = 10;
    simCtx->tiout = 10; simCtx->particleConsoleOutputFreq = simCtx->tiout;
    simCtx->StartTime = 0.0; simCtx->dt = 0.001;
    simCtx->OnlySetup = PETSC_FALSE;
    simCtx->continueMode = PETSC_FALSE;
    simCtx->logviewer = NULL;
    strcpy(simCtx->eulerianSource,"solve");
    strcpy(simCtx->restart_dir,"restart");
    strcpy(simCtx->output_dir,"output");
    strcpy(simCtx->log_dir,"logs");
    strcpy(simCtx->analysis_dir,"output/analysis/metrics");
    simCtx->_io_context_buffer[0] = '\0';
    simCtx->current_io_directory = NULL;
    simCtx->checkpointGeometrySHA256[0] = '\0';
    simCtx->checkpointGeometryHashReady = PETSC_FALSE;
    simCtx->restartHistoryAvailable = PETSC_FALSE;
 
    // --- Group 3: High-Level Physics & Model Selection Flags ---
    simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
    simCtx->invicid = 0;
    simCtx->TwoD = 0; simCtx->moveframe = 0; simCtx->rotateframe = 0;
    strcpy(simCtx->AnalyticalSolutionType,"TGV3D");
 
    // --- Group 4: Immersed-body flux corrections (refused; see Group 3 parsing) ---
    simCtx->MHV=0; simCtx->LV=0;
 
    // --- Group 5: Solver & Numerics Parameters ---
    simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK; simCtx->mom_max_pseudo_steps = 50;
    simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor = 1.1; // raised from 1.05; less aggressive rejection
    simCtx->mom_atol = 1e-7; simCtx->mom_rtol = 1e-4;
    /* Residual-based convergence is the default. Leaving both at 0.0 selected the
     * update-only branch, where |dU| <= mom_atol can pass purely because dtau collapsed
     * (|dU| ~ dtau*|R|), converging on a state that does not satisfy the equations.
     * Shipped configs all set these explicitly; these defaults protect a minimal
     * user-written config. Set both non-positive to opt back out deliberately. */
    simCtx->mom_resid_atol = 1e-8; simCtx->mom_resid_rtol = 1e-3;
    simCtx->imp_stol = 1.e-8;
    simCtx->mglevels = 3;
    simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
    simCtx->poisson = 0;
    simCtx->central=0;
    /* pseudo_cfl and its bounds are now dimensionless Courant numbers: CFL = dtau * lambda_max,
       where lambda_max is the global spectral radius computed at each physical timestep.
       Stable range for 4-stage Jameson RK: ~0–2.83. Initial 0.5 gives a comfortable margin. */
    simCtx->ren = 100.0; simCtx->pseudo_cfl = 0.5;
    simCtx->max_pseudo_cfl = 2.0; simCtx->min_pseudo_cfl = 0.001;
    simCtx->pseudo_cfl_reduction_factor = 0.75;
    simCtx->pseudo_cfl_growth_factor = 1.1; // raised from 1.0; controller can now increase CFL
    simCtx->no_pseudo_cfl_backtrack = PETSC_FALSE;
    simCtx->mom_ratio_ema_alpha = 0.3;   /* moderate smoothing; set to 1.0 to recover original raw-ratio behavior */
    simCtx->mom_last_converged  = PETSC_TRUE;
    simCtx->mom_last_lambda_max = 0.0;   /* populated after first momentum solve */
    simCtx->mom_nk_monitor_history = PETSC_FALSE;
    simCtx->ps_ksp_pic_monitor_true_residual = PETSC_FALSE;
    simCtx->cdisx = 0.0; simCtx->cdisy = 0.0; simCtx->cdisz = 0.0;
    simCtx->initialConditionMode = IC_MODE_ZERO;
    simCtx->initialConditionField = IC_FIELD_UCAT;
    strcpy(simCtx->initialConditionDirectory, "config/initial_condition");
    simCtx->InitialConstantContra.x = 0.0;
    simCtx->InitialConstantContra.y = 0.0;
    simCtx->InitialConstantContra.z = 0.0;
    simCtx->flowDirection = FLOW_DIR_UNSET;
    simCtx->icVelocityPhysical = 0.0;
    simCtx->AnalyticalUniformVelocity.x = 0.0;
    simCtx->AnalyticalUniformVelocity.y = 0.0;
    simCtx->AnalyticalUniformVelocity.z = 0.0;
    simCtx->fieldStatisticsEnabled = PETSC_FALSE;
    simCtx->fieldStatisticsWindowCount = 0;
    simCtx->fieldStatisticsWindows = NULL;
    simCtx->statisticsConsoleOutputFreq = 0;
    simCtx->fieldStatisticsContinue = PETSC_FALSE;
    simCtx->solutionConvergenceEnabled = PETSC_TRUE;
    simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
    simCtx->solutionConvergencePeriodSteps = 0;
    simCtx->solutionConvergenceWindowSteps = 0;
    simCtx->solutionConvergenceSamplesRecorded = 0;
    simCtx->solutionConvergenceMeanSpeedHistory = NULL;
    simCtx->solutionConvergenceMeanKEHistory = NULL;
    simCtx->verificationDiffusivity.enabled = PETSC_FALSE;
    strcpy(simCtx->verificationDiffusivity.mode, "");
    strcpy(simCtx->verificationDiffusivity.profile, "");
    simCtx->verificationDiffusivity.gamma0 = 0.0;
    simCtx->verificationDiffusivity.slope_x = 0.0;
 
    // --- Group 6: Physical & Geometric Parameters ---
    simCtx->Flux_in = 1.0; simCtx->angle = 0.0;
    simCtx->max_angle = -54. * 3.1415926 / 180.;
    simCtx->CMx_c=0.0; simCtx->CMy_c=0.0; simCtx->CMz_c=0.0;
    simCtx->wall_roughness_height = 1e-16;
    simCtx->schmidt_number = 1.0; simCtx->Turbulent_schmidt_number = 0.7;
    simCtx->iem_constant = 2.0;
 
    // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
    simCtx->block_number = 1; simCtx->inletprofile = 1;
    simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
    simCtx->pseudo_periodic = 0;
    strcpy(simCtx->grid_file, "config/grid.run");
    simCtx->generate_grid = PETSC_FALSE;
    simCtx->da_procs_x = PETSC_DECIDE;
    simCtx->da_procs_y = PETSC_DECIDE;
    simCtx->da_procs_z = PETSC_DECIDE;
    simCtx->num_bcs_files = 1;
    ierr = PetscMalloc1(1, &simCtx->bcs_files); CHKERRQ(ierr);
    ierr = PetscStrallocpy("config/bcs.run", &simCtx->bcs_files[0]); CHKERRQ(ierr);
    simCtx->FluxInSum = 0.0; simCtx->FluxOutSum = 0.0; simCtx->Fluxsum = 0.0;
    simCtx->drivingForceMagnitude = 0.0, simCtx->forceScalingFactor = 1.8;
    simCtx->drivingForceStep = -1;
    simCtx->drivenFluxMeasured = 0.0; simCtx->drivenFluxArea = 0.0;
    simCtx->targetVolumetricFlux  = 0.0;
    simCtx->drivenFluxTargetLatched = PETSC_FALSE;
    simCtx->bulkVelocityCorrection = 0.0;
    simCtx->boundaryVelocityCorrection = 0.0;
    simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
    simCtx->ccc = 0;
    simCtx->ratio = 0.0;
    
    
    // --- Group 8: Turbulence Modeling (LES) ---
    simCtx->les = NO_LES_MODEL;
    simCtx->wallfunction = 0; simCtx->les_gradient_model = 0;
    ierr = LESConfigSetDefaults(&simCtx->les_config); CHKERRQ(ierr);
 
    // --- Group 9: Particle / DMSwarm Data & Settings ---
    simCtx->np = 0;
    simCtx->particleRandomSeed = 12345;
    simCtx->dm_swarm = NULL; simCtx->bboxlist = NULL;
    simCtx->ParticleInitialization = PARTICLE_INIT_SURFACE_RANDOM;
    simCtx->interpolationMethod = INTERP_TRILINEAR;
    strcpy(simCtx->particleRestartMode,"load");
    simCtx->particlesLostLastStep = 0;
    simCtx->particlesLostCumulative = 0;
    simCtx->particlesLostScalarLastStep = 0.0;
    simCtx->particlesMigratedLastStep = 0;
    simCtx->occupiedCellCount = 0;
    simCtx->particleLoadImbalance = 0.0;
    simCtx->migrationPassesLastStep = 0;
    simCtx->searchMetrics.searchAttempts = 0;
    simCtx->searchMetrics.searchPopulation = 0;
    simCtx->searchMetrics.searchLocatedCount = 0;
    simCtx->searchMetrics.searchLostCount = 0;
    simCtx->searchMetrics.traversalStepsSum = 0;
    simCtx->searchMetrics.reSearchCount = 0;
    simCtx->searchMetrics.maxTraversalSteps = 0;
    simCtx->searchMetrics.maxTraversalFailCount = 0;
    simCtx->searchMetrics.tieBreakCount = 0;
    simCtx->searchMetrics.boundaryClampCount = 0;
    simCtx->searchMetrics.bboxGuessSuccessCount = 0;
    simCtx->searchMetrics.bboxGuessFallbackCount = 0;
    simCtx->searchMetrics.maxParticlePassDepth = 0;
    simCtx->searchMetrics.currentSettlementPass = 0;
    simCtx->BrownianMotionRNG = NULL;
    simCtx->C_IEM = 2.0;
 
    // --- Group 10: Immersed Boundary & FSI Data Object Pointers ---
    simCtx->ibm = NULL; simCtx->ibmv = NULL; simCtx->fsi = NULL;
 
    // --- Group 11: Logging and Custom Configuration ---
    strcpy(simCtx->allowedFile, "config/whitelist.run");
    simCtx->useCfg = PETSC_FALSE;
    simCtx->allowedFuncs = NULL;
    simCtx->nAllowed = 0;
    simCtx->LoggingFrequency = 10;
    simCtx->poissonSourceImbalance = 0.0;
    simCtx->MaxDiv = 0.0;
    simCtx->MaxDivFlatArg = 0; simCtx->MaxDivx = 0; simCtx->MaxDivy = 0; simCtx->MaxDivz = 0;
    strcpy(simCtx->profilingSelectedFuncsFile, "config/profile.run");
    simCtx->useProfilingSelectedFuncsCfg =  PETSC_FALSE;
    simCtx->profilingSelectedFuncs = NULL;
    simCtx->nProfilingSelectedFuncs = 0;
    strcpy(simCtx->profilingTimestepMode, "selected");
    strcpy(simCtx->profilingTimestepFile, "Profiling_Timestep_Summary.csv");
    simCtx->profilingFinalSummary = PETSC_TRUE;
    simCtx->walltimeGuardEnabled = PETSC_FALSE;
    simCtx->walltimeGuardActive = PETSC_FALSE;
    simCtx->walltimeGuardWarmupSteps = 10;
    simCtx->walltimeGuardMultiplier = 2.0;
    simCtx->walltimeGuardMinSeconds = 60.0;
    simCtx->walltimeGuardEstimatorAlpha = 0.35;
    simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
    simCtx->walltimeGuardLimitSeconds = 0.0;
    simCtx->walltimeGuardCompletedSteps = 0;
    simCtx->walltimeGuardWarmupTotalSeconds = 0.0;
    simCtx->walltimeGuardWarmupAverageSeconds = 0.0;
    simCtx->walltimeGuardHasEWMA = PETSC_FALSE;
    simCtx->walltimeGuardEWMASeconds = 0.0;
    simCtx->walltimeGuardLatestStepSeconds = 0.0;
    simCtx->runtimeMemoryLogEnabled = PETSC_TRUE;
    strcpy(simCtx->runtimeMemoryLogFile, "Runtime_Memory.log");
    simCtx->runtimeMemoryLogStarted = PETSC_FALSE;
    simCtx->runtimeMemoryLogHasPrevious = PETSC_FALSE;
    simCtx->runtimeMemoryLogPreviousProcessMB = 0.0;
    // --- Group 11: Post-Processing Information ---
    strcpy(simCtx->PostprocessingControlFile, "config/post.run");
    ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
 
    // === 2. Get MPI Info and Handle Config File =============================
    // -- Group 1:  Parallelism & MPI Information
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &simCtx->rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &simCtx->size); CHKERRQ(ierr);
    
 // First, check if the -control_file argument was provided by the user/script.
    ierr = PetscOptionsGetString(NULL, NULL, "-control_file", control_filename, sizeof(control_filename), &control_flg); CHKERRQ(ierr);
 
    // If the flag is NOT present or the filename is empty, abort with a helpful error.
    if (!control_flg || strlen(control_filename) == 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "\n\n*** MANDATORY ARGUMENT MISSING ***\n"
                "The -control_file argument was not provided.\n"
                "This program must be launched with a configuration file.\n"
                "Example: mpiexec -n 4 ./simulator -control_file /path/to/your/config.control\n"
                "This is typically handled automatically by the 'picurv' script.\n");
    }
 
    // At this point, we have a valid filename. Attempt to load it.
    LOG(GLOBAL, LOG_INFO, "Loading mandatory configuration from: %s\n", control_filename);
    ierr = PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, control_filename, PETSC_FALSE);
    if (ierr == PETSC_ERR_FILE_OPEN) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "The specified control file was not found or could not be opened: %s", control_filename);
    }
    CHKERRQ(ierr);
 
    {
        PetscBool legacy_averaging = PETSC_FALSE;
        ierr = PetscOptionsHasName(NULL, NULL, "-averaging", &legacy_averaging); CHKERRQ(ierr);
        PetscCheck(!legacy_averaging, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                   "Legacy -averaging was removed. Use instantaneous output and offline "
                   "postprocessing until the replacement field-statistics pipeline is available.");
    }
 
    // === 3. A Configure Logging System ========================================
    // This logic determines the logging configuration and STORES it in simCtx for
    // later reference and cleanup.
    ierr = PetscOptionsGetString(NULL, NULL, "-whitelist_config_file", simCtx->allowedFile, PETSC_MAX_PATH_LEN, &simCtx->useCfg); CHKERRQ(ierr);
    
    if (simCtx->useCfg) {
        ierr = LoadAllowedFunctionsFromFile(simCtx->allowedFile, &simCtx->allowedFuncs, &simCtx->nAllowed);
        if (ierr) {
            // Use direct PetscPrintf as logging system isn't fully active yet.
            PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load allowed functions from '%s'. Falling back to default list.\n", __func__, simCtx->allowedFile);
            simCtx->useCfg = PETSC_FALSE; // Mark as failed.
            ierr = 0; // Clear the error to allow fallback.
        } else if (simCtx->nAllowed == 0) {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                    "Whitelist config file '%s' is empty. Omit -whitelist_config_file to use the default allow-list, or list at least one function.",
                    simCtx->allowedFile);
        }
    }
    if (!simCtx->useCfg) {
        // Fallback to default logging functions if no file was used or if loading failed.
        simCtx->nAllowed = 2;
        ierr = PetscMalloc1(simCtx->nAllowed, &simCtx->allowedFuncs); CHKERRQ(ierr);
        ierr = PetscStrallocpy("main", &simCtx->allowedFuncs[0]); CHKERRQ(ierr);
        ierr = PetscStrallocpy("CreateSimulationContext", &simCtx->allowedFuncs[1]); CHKERRQ(ierr);
    }
    
    // Activate the configuration by passing it to the logging module's setup function.
    set_allowed_functions((const char**)simCtx->allowedFuncs, (size_t)simCtx->nAllowed);
    
    // Now that the logger is configured, we can use it.
    LOG_ALLOW_SYNC(LOCAL, LOG_INFO, "Context created. Initializing on rank %d of %d.\n", simCtx->rank, simCtx->size);
    print_log_level(); // This will now correctly reflect the LOG_LEVEL environment variable.
   
    // === 3.B Configure Profiling System ========================================
    ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_mode", simCtx->profilingTimestepMode, sizeof(simCtx->profilingTimestepMode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_file", simCtx->profilingTimestepFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-profiling_final_summary", &simCtx->profilingFinalSummary, NULL); CHKERRQ(ierr);
    if (strcmp(simCtx->profilingTimestepMode, "off") != 0 &&
        strcmp(simCtx->profilingTimestepMode, "selected") != 0 &&
        strcmp(simCtx->profilingTimestepMode, "all") != 0) {
        PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Unknown profiling timestep mode '%s'. Falling back to 'selected'.\n", __func__, simCtx->profilingTimestepMode);
        strcpy(simCtx->profilingTimestepMode, "selected");
    }
 
    if (strcmp(simCtx->profilingTimestepMode, "selected") == 0) {
        ierr = PetscOptionsGetString(NULL, NULL, "-profile_config_file", simCtx->profilingSelectedFuncsFile, PETSC_MAX_PATH_LEN, &simCtx->useProfilingSelectedFuncsCfg); CHKERRQ(ierr);
        if (simCtx->useProfilingSelectedFuncsCfg) {
            ierr = LoadAllowedFunctionsFromFile(simCtx->profilingSelectedFuncsFile, &simCtx->profilingSelectedFuncs, &simCtx->nProfilingSelectedFuncs);
            if (ierr) {
                PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load selected profiling functions from '%s'. Falling back to default list.\n", __func__, simCtx->profilingSelectedFuncsFile);
                simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
                ierr = 0;
            }
        }
        if (!simCtx->useProfilingSelectedFuncsCfg) {
            // Fallback to a hardcoded default list if no file was provided or loading failed.
            simCtx->nProfilingSelectedFuncs = 4;
            ierr = PetscMalloc1(simCtx->nProfilingSelectedFuncs, &simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
            ierr = PetscStrallocpy("FlowSolver", &simCtx->profilingSelectedFuncs[0]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("AdvanceSimulation", &simCtx->profilingSelectedFuncs[1]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("LocateAllParticlesInGrid", &simCtx->profilingSelectedFuncs[2]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("InterpolateAllFieldsToSwarm", &simCtx->profilingSelectedFuncs[3]); CHKERRQ(ierr);
        }
    }
 
    // Initialize the profiling system with the current updated simulation context.
    ierr = ProfilingInitialize(simCtx); CHKERRQ(ierr);
 
    // === 4. Parse All Command Line Options ==================================
    LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing command-line options...\n");
 
    // --- Group 2
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 2: Simulation Control,Time and I/O.\n");
    // Read the physical time to start from.
    // The default is already 0.0, so this will only be non-zero if the user provides it.
    ierr = PetscOptionsGetInt(NULL, NULL, "-start_step", &simCtx->StartStep, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,NULL, "-totalsteps", &simCtx->StepsToRun, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-only_setup", &simCtx->OnlySetup, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-continue_mode", &simCtx->continueMode, NULL); CHKERRQ(ierr);
    /* Resuming in the same run directory continues accumulated statistics by
     * default; branching with --restart-from must opt in, because the branch may
     * follow a different physical trajectory than the samples already collected. */
    simCtx->fieldStatisticsContinue = simCtx->continueMode;
    ierr = PetscOptionsGetBool(NULL, NULL, "-field_statistics_continue",
                               &simCtx->fieldStatisticsContinue, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-dt", &simCtx->dt, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-tio", &simCtx->tiout, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-particle_console_output_freq", &simCtx->particleConsoleOutputFreq, &particle_console_output_freq_flg); CHKERRQ(ierr);
    if (!particle_console_output_freq_flg) {
        simCtx->particleConsoleOutputFreq = simCtx->tiout;
    }
    ierr = PetscOptionsGetString(NULL,NULL,"-euler_field_source",simCtx->eulerianSource,sizeof(simCtx->eulerianSource),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-output_dir",simCtx->output_dir,sizeof(simCtx->output_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-restart_dir",simCtx->restart_dir,sizeof(simCtx->restart_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-log_dir",simCtx->log_dir,sizeof(simCtx->log_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-analysis_dir",simCtx->analysis_dir,sizeof(simCtx->analysis_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-walltime_guard_enabled", &simCtx->walltimeGuardEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-walltime_guard_warmup_steps", &simCtx->walltimeGuardWarmupSteps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_multiplier", &simCtx->walltimeGuardMultiplier, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-runtime_memory_log_enabled", &simCtx->runtimeMemoryLogEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-runtime_memory_log_file", simCtx->runtimeMemoryLogFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_min_seconds", &simCtx->walltimeGuardMinSeconds, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_estimator_alpha", &simCtx->walltimeGuardEstimatorAlpha, NULL); CHKERRQ(ierr);
 
    if (simCtx->walltimeGuardWarmupSteps <= 0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_warmup_steps: %d. Must be > 0.", simCtx->walltimeGuardWarmupSteps);
    }
    if (simCtx->walltimeGuardMultiplier <= 0.0 || simCtx->walltimeGuardMultiplier > 5.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_multiplier: %.6f. Must be in (0, 5].", (double)simCtx->walltimeGuardMultiplier);
    }
    if (simCtx->walltimeGuardMinSeconds <= 0.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_min_seconds: %.6f. Must be > 0.", (double)simCtx->walltimeGuardMinSeconds);
    }
    if (simCtx->walltimeGuardEstimatorAlpha <= 0.0 || simCtx->walltimeGuardEstimatorAlpha > 1.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_estimator_alpha: %.6f. Must be in (0, 1].", (double)simCtx->walltimeGuardEstimatorAlpha);
    }
 
    if(strcmp(simCtx->eulerianSource,"solve")!= 0 && strcmp(simCtx->eulerianSource,"load") != 0 && strcmp(simCtx->eulerianSource,"analytical")!=0){
      SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Invalid value for -euler_field_source. Must be 'load','analytical' or 'solve'. You provided '%s'.",simCtx->eulerianSource);
    }
    if (simCtx->walltimeGuardEnabled) {
      const char *job_start_env = getenv("PICURV_JOB_START_EPOCH");
      const char *limit_env = getenv("PICURV_WALLTIME_LIMIT_SECONDS");
      PetscBool job_start_ok = RuntimeWalltimeGuardParsePositiveSeconds(job_start_env, &simCtx->walltimeGuardJobStartEpochSeconds);
      PetscBool limit_ok = RuntimeWalltimeGuardParsePositiveSeconds(limit_env, &simCtx->walltimeGuardLimitSeconds);
 
      if (!job_start_ok || !limit_ok) {
        simCtx->walltimeGuardActive = PETSC_FALSE;
        simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
        simCtx->walltimeGuardLimitSeconds = 0.0;
        LOG_ALLOW(
          GLOBAL,
          LOG_WARNING,
          "Runtime walltime guard enabled but %s/%s are missing or invalid. Falling back to external shutdown signals only.\n",
          "PICURV_JOB_START_EPOCH",
          "PICURV_WALLTIME_LIMIT_SECONDS"
        );
      } else {
        simCtx->walltimeGuardActive = PETSC_TRUE;
      }
    }
 
    //  --- Group 3
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 3: High-Level Physics & Model Selection Flags\n");    
    ierr = PetscOptionsGetInt(NULL, NULL, "-imm", &simCtx->immersed, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-fsi", &simCtx->movefsi, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rfsi", &simCtx->rotatefsi, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-inv", &simCtx->invicid, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-TwoD", &simCtx->TwoD, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-analytical_type",simCtx->AnalyticalSolutionType,sizeof(simCtx->AnalyticalSolutionType),NULL);CHKERRQ(ierr);
 
    /* Immersed boundaries, moving bodies and moving reference frames are planned, not
       implemented. No body geometry is ever loaded, and the IBM interpolation the
       momentum solvers would call is commented out and defined nowhere, so -imm only
       reconfigures the Poisson solve around a body that is not there. The moving-frame
       convection branch in ComputeRHS is commented out as well, so -mframe/-rframe would
       drop the convective term. Each would run a different problem without saying so;
       the conductor refuses the YAML switches, and this refuses the same flags arriving
       through a PETSc passthrough. */
    PetscCheck(!simCtx->immersed && !simCtx->movefsi && !simCtx->rotatefsi, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Immersed boundaries and moving bodies (-imm, -fsi, -rfsi) are planned, not implemented.");
    PetscCheck(!simCtx->moveframe && !simCtx->rotateframe, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Moving and rotating reference frames (-mframe, -rframe) are planned, not implemented: "
               "their convection branch is absent, so the convective term would be dropped.");
 
    //  --- Group 4
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Immersed-Body Flux Corrections and Driven-Flow Controls\n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
    /* The heart-valve volume-flux corrections balance flux across an immersed body, which
       cannot exist while immersed boundaries are unimplemented (see the Group 3 check). */
    PetscCheck(!simCtx->MHV && !simCtx->LV, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "The -mhv/-lv immersed-body flux corrections need immersed boundaries, which are planned, not implemented.");
    ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_initial_force",&simCtx->drivingForceMagnitude,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_scaling_factor",&simCtx->forceScalingFactor,NULL);CHKERRQ(ierr);
    //  --- Group 5
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 5: Solver & Numerics Parameters \n");
    char mom_solver_type_char[PETSC_MAX_PATH_LEN];
    char solution_convergence_mode_char[PETSC_MAX_PATH_LEN];
    PetscBool mom_solver_type_flg = PETSC_FALSE;
    PetscBool solution_convergence_mode_flg = PETSC_FALSE;
    ierr = PetscOptionsGetString(NULL, NULL, "-mom_solver_type", mom_solver_type_char, sizeof(mom_solver_type_char), &mom_solver_type_flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mom_max_pseudo_steps", &simCtx->mom_max_pseudo_steps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_atol", &simCtx->mom_atol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_rtol", &simCtx->mom_rtol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_atol", &simCtx->mom_resid_atol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_rtol", &simCtx->mom_resid_rtol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-imp_stol", &simCtx->imp_stol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-central", &simCtx->central, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-solution_convergence_mode",
                                 solution_convergence_mode_char, sizeof(solution_convergence_mode_char),
                                 &solution_convergence_mode_flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-solution_convergence_enabled", &simCtx->solutionConvergenceEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_period_steps", &simCtx->solutionConvergencePeriodSteps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_window_steps", &simCtx->solutionConvergenceWindowSteps, NULL); CHKERRQ(ierr);
    /* Field statistics resolve here rather than inline, because a variable-arity
     * window list needs its own parse site; see statistics_config.c. It must
     * precede CreateAndInitializeAllVectors, which sizes the per-window
     * accumulators from the window count resolved here. */
    ierr = ParseFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
 
    // Keep parser acceptance aligned with the enum and FlowSolver dispatch.
    if (mom_solver_type_flg) {
        if(strcmp(mom_solver_type_char, "DUALTIME_PICARD_JAMESON_RK") == 0 ||
           strcmp(mom_solver_type_char, "DUALTIME_PICARD_RK4") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK;
        } else if (strcmp(mom_solver_type_char, "EXPLICIT_RK") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_EXPLICIT_RK;
        } else if (strcmp(mom_solver_type_char, "newton_krylov") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_NEWTON_KRYLOV;
        } else {
            LOG(GLOBAL, LOG_ERROR, "Invalid value for -mom_solver_type: '%s'. Valid options are: 'DUALTIME_PICARD_JAMESON_RK', 'EXPLICIT_RK', 'newton_krylov'.\n", mom_solver_type_char);
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -mom_solver_type: '%s'.", mom_solver_type_char);
        }
    }
 
    if (solution_convergence_mode_flg) {
        if (strcmp(solution_convergence_mode_char, "STEADY_DETERMINISTIC") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
        } else if (strcmp(solution_convergence_mode_char, "PERIODIC_DETERMINISTIC") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC;
        } else if (strcmp(solution_convergence_mode_char, "STATISTICAL_STEADY") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STATISTICAL_STEADY;
        } else if (strcmp(solution_convergence_mode_char, "TRANSIENT") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_TRANSIENT;
        } else {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                    "Invalid value for -solution_convergence_mode: '%s'.", solution_convergence_mode_char);
        }
    }
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
        simCtx->solutionConvergencePeriodSteps <= 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "solution convergence mode PERIODIC_DETERMINISTIC requires -solution_convergence_period_steps > 0.");
    }
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY &&
        simCtx->solutionConvergenceWindowSteps <= 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "solution convergence mode STATISTICAL_STEADY requires -solution_convergence_window_steps > 0.");
    }
 
    // --- Multigrid Options ---
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_level", &simCtx->mglevels, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_pre_it", &simCtx->mg_preItr, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_post_it", &simCtx->mg_poItr, NULL); CHKERRQ(ierr);
 
    // --- Other Solver Options ---
    ierr = PetscOptionsGetInt(NULL, NULL, "-poisson", &simCtx->poisson, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ren", &simCtx->ren, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl", &simCtx->pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-max_pseudo_cfl", &simCtx->max_pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-min_pseudo_cfl", &simCtx->min_pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_reduction_factor", &simCtx->pseudo_cfl_reduction_factor, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_growth_factor", &simCtx->pseudo_cfl_growth_factor, NULL); CHKERRQ(ierr);
    // Read the deprecated RK4 spelling first so the canonical Jameson option wins if both are present.
    ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_rk4_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_jameson_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-no_pseudo_cfl_backtrack", &simCtx->no_pseudo_cfl_backtrack, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_ratio_ema_alpha", &simCtx->mom_ratio_ema_alpha, NULL); CHKERRQ(ierr);
    if (simCtx->min_pseudo_cfl <= 0.0 ||
        simCtx->pseudo_cfl < simCtx->min_pseudo_cfl ||
        simCtx->pseudo_cfl > simCtx->max_pseudo_cfl) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "Pseudo-CFL controls require 0 < minimum <= initial <= maximum.");
    }
    if (simCtx->pseudo_cfl_growth_factor < 1.0 ||
        simCtx->pseudo_cfl_reduction_factor <= 0.0 ||
        simCtx->pseudo_cfl_reduction_factor >= 1.0 ||
        simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor < 1.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "Pseudo-CFL controls require growth_factor >= 1, 0 < reduction_factor < 1, and noise allowance >= 1.");
    }
    if (simCtx->mom_ratio_ema_alpha < 0.0 || simCtx->mom_ratio_ema_alpha > 1.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "-mom_ratio_ema_alpha must be in [0, 1].");
    }
    ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &simCtx->ps_ksp_pic_monitor_true_residual); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-mom_nk_pic_monitor", &simCtx->mom_nk_monitor_history, NULL); CHKERRQ(ierr);
    {
        PetscInt ic_mode = (PetscInt)simCtx->initialConditionMode;
        PetscInt ic_field = (PetscInt)simCtx->initialConditionField;
        ierr = PetscOptionsGetInt(NULL, NULL, "-finit", &ic_mode, NULL); CHKERRQ(ierr);
        ierr = PetscOptionsGetInt(NULL, NULL, "-ic_field", &ic_field, NULL); CHKERRQ(ierr);
        simCtx->initialConditionMode = (InitialConditionMode)ic_mode;
        simCtx->initialConditionField = (InitialConditionField)ic_field;
    }
    ierr = PetscOptionsGetString(NULL, NULL, "-ic_dir", simCtx->initialConditionDirectory,
                                 sizeof(simCtx->initialConditionDirectory), NULL); CHKERRQ(ierr);
    if (simCtx->initialConditionMode < IC_MODE_ZERO || simCtx->initialConditionMode > IC_MODE_FILE) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
              "Invalid value for -finit. Expected an initial-condition mode in [0,4], got %d.",
              simCtx->initialConditionMode);
    }
    if (simCtx->initialConditionField < IC_FIELD_UCAT || simCtx->initialConditionField > IC_FIELD_UCONT) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
              "Invalid value for -ic_field. Expected 0 (Ucat) or 1 (Ucont), got %d.",
              simCtx->initialConditionField);
    }
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_x", &simCtx->InitialConstantContra.x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_y", &simCtx->InitialConstantContra.y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_z", &simCtx->InitialConstantContra.z, NULL); CHKERRQ(ierr);
    {
        PetscInt  fd_int = (PetscInt)FLOW_DIR_UNSET;
        PetscBool fd_set = PETSC_FALSE;
        ierr = PetscOptionsGetInt(NULL, NULL, "-flow_direction", &fd_int, &fd_set); CHKERRQ(ierr);
        if (fd_set) simCtx->flowDirection = (FlowDirection)fd_int;
    }
    ierr = PetscOptionsGetReal(NULL, NULL, "-ic_velocity_physical", &simCtx->icVelocityPhysical, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_u", &simCtx->AnalyticalUniformVelocity.x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_v", &simCtx->AnalyticalUniformVelocity.y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_w", &simCtx->AnalyticalUniformVelocity.z, NULL); CHKERRQ(ierr);
    PetscBool verification_scalar_value_set = PETSC_FALSE;
    PetscBool verification_scalar_phi0_set = PETSC_FALSE;
    PetscBool verification_scalar_slope_x_set = PETSC_FALSE;
    PetscBool verification_scalar_amplitude_set = PETSC_FALSE;
    PetscBool verification_scalar_kx_set = PETSC_FALSE;
    PetscBool verification_scalar_ky_set = PETSC_FALSE;
    PetscBool verification_scalar_kz_set = PETSC_FALSE;
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_mode",
                                 simCtx->verificationDiffusivity.mode,
                                 sizeof(simCtx->verificationDiffusivity.mode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_profile",
                                 simCtx->verificationDiffusivity.profile,
                                 sizeof(simCtx->verificationDiffusivity.profile), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_gamma0",
                               &simCtx->verificationDiffusivity.gamma0, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_slope_x",
                               &simCtx->verificationDiffusivity.slope_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_mode",
                                 simCtx->verificationScalar.mode,
                                 sizeof(simCtx->verificationScalar.mode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_profile",
                                 simCtx->verificationScalar.profile,
                                 sizeof(simCtx->verificationScalar.profile), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_value",
                               &simCtx->verificationScalar.value, &verification_scalar_value_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_phi0",
                               &simCtx->verificationScalar.phi0, &verification_scalar_phi0_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_slope_x",
                               &simCtx->verificationScalar.slope_x, &verification_scalar_slope_x_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_amplitude",
                               &simCtx->verificationScalar.amplitude, &verification_scalar_amplitude_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kx",
                               &simCtx->verificationScalar.kx, &verification_scalar_kx_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_ky",
                               &simCtx->verificationScalar.ky, &verification_scalar_ky_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kz",
                               &simCtx->verificationScalar.kz, &verification_scalar_kz_set); CHKERRQ(ierr);
    simCtx->verificationDiffusivity.enabled =
      (PetscBool)(simCtx->verificationDiffusivity.mode[0] != '\0' ||
                  simCtx->verificationDiffusivity.profile[0] != '\0');
    simCtx->verificationScalar.enabled =
      (PetscBool)(simCtx->verificationScalar.mode[0] != '\0' ||
                  simCtx->verificationScalar.profile[0] != '\0');
    if (simCtx->verificationDiffusivity.enabled) {
      if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                "verification diffusivity overrides require -euler_field_source \"analytical\".");
      }
      if (strcmp(simCtx->verificationDiffusivity.mode, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_diffusivity_mode '%s'. Only 'analytical' is supported.",
                simCtx->verificationDiffusivity.mode);
      }
      if (strcmp(simCtx->verificationDiffusivity.profile, "LINEAR_X") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_diffusivity_profile '%s'. Only 'LINEAR_X' is supported.",
                simCtx->verificationDiffusivity.profile);
      }
    }
    if (simCtx->verificationScalar.enabled) {
      if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                "verification scalar overrides require -euler_field_source \"analytical\".");
      }
      if (strcmp(simCtx->verificationScalar.mode, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_scalar_mode '%s'. Only 'analytical' is supported.",
                simCtx->verificationScalar.mode);
      }
      if (strcmp(simCtx->verificationScalar.profile, "CONSTANT") == 0) {
        if (!verification_scalar_value_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile CONSTANT requires -verification_scalar_value.");
        }
      } else if (strcmp(simCtx->verificationScalar.profile, "LINEAR_X") == 0) {
        if (!verification_scalar_phi0_set || !verification_scalar_slope_x_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile LINEAR_X requires -verification_scalar_phi0 and -verification_scalar_slope_x.");
        }
      } else if (strcmp(simCtx->verificationScalar.profile, "SIN_PRODUCT") == 0) {
        if (!verification_scalar_amplitude_set || !verification_scalar_kx_set ||
            !verification_scalar_ky_set || !verification_scalar_kz_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile SIN_PRODUCT requires -verification_scalar_amplitude, -verification_scalar_kx, -verification_scalar_ky, and -verification_scalar_kz.");
        }
      } else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_scalar_profile '%s'. Supported profiles: CONSTANT, LINEAR_X, SIN_PRODUCT.",
                simCtx->verificationScalar.profile);
      }
    }
    // NOTE: cdisx,cdisy,cdisz haven't been parsed, add if necessary.
 
     //  --- Group 6
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 6: Physical & Geometric Parameters \n");   
    ierr = PetscOptionsGetReal(NULL,NULL,"-schmidt_number",&simCtx->schmidt_number,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-turb_schmidt_number",&simCtx->Turbulent_schmidt_number,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-iem_constant",&simCtx->iem_constant,NULL);CHKERRQ(ierr);
    /* Zero is valid: it switches micromixing off, so a particle scalar is carried unchanged. */
    PetscCheck(simCtx->iem_constant >= 0.0 && !PetscIsInfOrNanReal(simCtx->iem_constant),
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-iem_constant must be a non-negative finite number (got %g).", (double)simCtx->iem_constant);
    ierr = PetscOptionsGetReal(NULL,NULL,"-wall_roughness",&simCtx->wall_roughness_height,NULL);CHKERRQ(ierr);
    // NOTE: angle is not parsed in the original code, it set programmatically. We will follow that.
    // NOTE: max_angle is calculated based on other flags (like MHV) in the legacy code.
    // We will defer that logic to a later setup stage and not parse them directly.
    // The Scaling Information is calculated here
    ierr = ParseScalingInformation(simCtx); CHKERRQ(ierr);
 
     //  --- Group 7
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 7: Grid, Domain, and Boundary Condition Settings \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-nblk", &simCtx->block_number, NULL); CHKERRQ(ierr); // This is also a modern option
    /* Every inter-block exchange (Block_Interface_U) is commented out in the momentum
       solvers and the initial condition, so blocks would be solved as isolated domains.
       Multi-block coupling is planned, not implemented. */
    PetscCheck(simCtx->block_number == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Multi-block domains are planned, not implemented: blocks are never coupled (got -nblk %" PetscInt_FMT ").",
               simCtx->block_number);
    ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
    // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
    ierr = PetscOptionsGetBool(NULL, NULL, "-grid", &simCtx->generate_grid, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-grid_file", simCtx->grid_file, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_x", &simCtx->da_procs_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_y", &simCtx->da_procs_y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_z", &simCtx->da_procs_z, NULL); CHKERRQ(ierr);
    // NOTE: pseudo_periodic was not parsed. We will omit its parsing call. 
    PetscBool bcs_flg;
    char      file_list_str[PETSC_MAX_PATH_LEN * 10]; // Buffer for comma-separated list
 
    ierr = PetscOptionsGetString(NULL, NULL, "-bcs_files", file_list_str, sizeof(file_list_str), &bcs_flg); CHKERRQ(ierr);
 
    if (bcs_flg) {
      LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found -bcs_files option, overriding default.\n");
    
      // A. Clean up the default memory we allocated in Phase 1.
      ierr = PetscFree(simCtx->bcs_files[0]); CHKERRQ(ierr);
      ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
      simCtx->num_bcs_files = 0;
      simCtx->bcs_files = NULL;
 
      // B. Parse the user-provided comma-separated list.
      char *token;
      char *str_copy;
      ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
 
      // First pass: count the number of files.
      token = strtok(str_copy, ",");
      while (token) {
        simCtx->num_bcs_files++;
        token = strtok(NULL, ",");
      }
      ierr = PetscFree(str_copy); CHKERRQ(ierr);
    
      // Second pass: allocate memory and store the filenames.
      ierr = PetscMalloc1(simCtx->num_bcs_files, &simCtx->bcs_files); CHKERRQ(ierr);
      ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
      token = strtok(str_copy, ",");
      for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
        ierr = PetscStrallocpy(token, &simCtx->bcs_files[i]); CHKERRQ(ierr);
        token = strtok(NULL, ",");
      }
      ierr = PetscFree(str_copy); CHKERRQ(ierr);
    }
    
 
     //  --- Group 8
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 8: Turbulence Modeling (LES) \n");
    // Seeded from the default already in simCtx: PetscOptionsGetInt leaves the target
    // untouched when -les is absent, and an uninitialized value would become the model.
    PetscInt temp_les_model = (PetscInt)simCtx->les;
    ierr = PetscOptionsGetInt(NULL, NULL, "-les", &temp_les_model, NULL); CHKERRQ(ierr);
    PetscCheck(temp_les_model >= NO_LES_MODEL && temp_les_model <= WALE, PETSC_COMM_WORLD,
               PETSC_ERR_ARG_OUTOFRANGE,
               "-les must be 0 (none), 1 (constant_smagorinsky), 2 (dynamic_smagorinsky), "
               "3 (vreman), or 4 (wale); received %" PetscInt_FMT ".", temp_les_model);
    simCtx->les = (LESModelType)temp_les_model;
    {
        /* RANS is planned, not implemented: the k-omega transport update was never
           written and its fields were never allocated. The selector and its dead
           hooks were removed on 2026-09-22 (see src/guide.md); this refuses the flag
           arriving through a PETSc passthrough. */
        PetscInt requested_rans = 0;
        ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &requested_rans, NULL); CHKERRQ(ierr);
        PetscCheck(requested_rans == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP,
                   "RANS closures (-rans) are planned, not implemented: no transport equation is "
                   "solved for them. Use an LES closure, or run laminar.");
    }
    ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
    PetscCheck(simCtx->wallfunction >= WALL_FUNCTION_NONE && simCtx->wallfunction <= WALL_FUNCTION_CABOT,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-wallfunction must be 0 (none), 1 (log_law), 2 (werner), or 3 (cabot); received %" PetscInt_FMT ".",
               simCtx->wallfunction);
    ierr = PetscOptionsGetInt(NULL, NULL, "-les_gradient_model", &simCtx->les_gradient_model, NULL); CHKERRQ(ierr);
    ierr = ParseLESConfiguration(simCtx); CHKERRQ(ierr);
 
     //  --- Group 9
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 9:  Particle / DMSwarm Data & Settings \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-numParticles", &simCtx->np, NULL); CHKERRQ(ierr);
    PetscInt temp_pinit = (PetscInt)PARTICLE_INIT_SURFACE_RANDOM;
    ierr = PetscOptionsGetInt(NULL, NULL, "-pinit", &temp_pinit, NULL); CHKERRQ(ierr);
    simCtx->ParticleInitialization = (ParticleInitializationType)temp_pinit;
    PetscInt temp_interp = (PetscInt)INTERP_TRILINEAR;
    ierr = PetscOptionsGetInt(NULL, NULL, "-interpolation_method", &temp_interp, NULL); CHKERRQ(ierr);
    simCtx->interpolationMethod = (InterpolationMethod)temp_interp;
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Interpolation method: %s\n",
              simCtx->interpolationMethod == INTERP_TRILINEAR ? "Trilinear (direct cell-center)" : "CornerAveraged (legacy)");
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_x", &simCtx->psrc_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_y", &simCtx->psrc_y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_z", &simCtx->psrc_z, NULL); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Particle initialization mode: %s. Point source: (%.6f, %.6f, %.6f)\n",
              ParticleInitializationToString(simCtx->ParticleInitialization),
              simCtx->psrc_x, simCtx->psrc_y, simCtx->psrc_z);
    ierr = PetscOptionsGetString(NULL,NULL,"-particle_restart_mode",simCtx->particleRestartMode,sizeof(simCtx->particleRestartMode),NULL); CHKERRQ(ierr);
    // Validation for Particle Restart Mode
    if (strcmp(simCtx->particleRestartMode, "load") != 0 && strcmp(simCtx->particleRestartMode, "init") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -particle_restart_mode. Must be 'load' or 'init'. You provided '%s'.", simCtx->particleRestartMode);
    }
    ierr = PetscOptionsGetInt(NULL, NULL, "-particle_random_seed", &simCtx->particleRandomSeed, NULL); CHKERRQ(ierr);
    PetscCheck(simCtx->particleRandomSeed >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-particle_random_seed must be non-negative (got %" PetscInt_FMT ").", simCtx->particleRandomSeed);
    ierr = InitializeBrownianRNG(simCtx); CHKERRQ(ierr);
    ierr = ParticleFieldPlanCreate(&simCtx->particleFieldPlan); CHKERRQ(ierr);
    PetscCheck(!simCtx->particleFieldPlan || simCtx->np > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "Particle field initial values are configured but the run has no particles.");
    PetscCheck(!simCtx->particleFieldPlan || !VerificationScalarOverrideActive(simCtx), PETSC_COMM_WORLD,
               PETSC_ERR_ARG_WRONG,
               "Particle field initial values cannot be combined with the verification scalar source, "
               "which prescribes Psi at every step.");
    // --- Group 10
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
 
    // --- Group 11
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 11: Top-Level Managers & Custom Configuration \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-logfreq", &simCtx->LoggingFrequency, NULL); CHKERRQ(ierr);
 
    if (simCtx->num_bcs_files != simCtx->block_number) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Number of BC files (%d) does not match number of blocks (%d). Use -bcs_files \"file1.dat,file2.dat,...\".", simCtx->num_bcs_files, simCtx->block_number);
    }
    
    // --- Group 12
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 12: Post-Processing Information.\n");
    // This logic determines the Post Processing configuration and STORES it in simCtx for later reference and cleanup.
    ierr = PetscOptionsGetString(NULL,NULL,"-postprocessing_config_file",simCtx->PostprocessingControlFile,PETSC_MAX_PATH_LEN,NULL); CHKERRQ(ierr);
    /* Parse post settings for both solver and post-processor binaries using the single pre-allocated pps object. */
    ierr = ParsePostProcessingSettings(simCtx); CHKERRQ(ierr);
 
    // === 5. Dependent Parameter Calculations ================================
    // Some parameters depend on others, so we calculate them here.
    simCtx->StartTime = (PetscReal)simCtx->StartStep*simCtx->dt;
    simCtx->ti = simCtx->StartTime;
    simCtx->step = simCtx->StartStep;
 
    // === 5. Log Summary and Finalize Setup ==================================
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "-- Console Output Functions [Total : %d] : --\n", simCtx->nAllowed);
    for (PetscInt i = 0; i < simCtx->nAllowed; ++i) {
      LOG_ALLOW(GLOBAL, LOG_DEBUG, "   [%2d] «%s»\n", i, simCtx->allowedFuncs[i]);
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration complete. Key parameters:\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Run mode: %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Time steps: %d (from %d to %d)\n", simCtx->StepsToRun, simCtx->StartStep, simCtx->StartStep + simCtx->StepsToRun);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Time step size (dt): %g\n", simCtx->dt);
    if (simCtx->tiout > 0) {
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Field/restart output cadence: every %d step(s)\n", simCtx->tiout);
    } else {
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Field/restart output cadence: DISABLED\n");
    }
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Immersed Boundary: %s\n", simCtx->immersed ? "ENABLED" : "DISABLED");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particles: %d\n", simCtx->np);
    if (simCtx->np > 0) {
      if (simCtx->particleConsoleOutputFreq > 0) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console cadence: every %d step(s)\n", simCtx->particleConsoleOutputFreq);
      } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console cadence: DISABLED\n");
      }
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console row subsampling: every %d particle(s)\n", simCtx->LoggingFrequency);
    }
    if (simCtx->StartStep > 0 && simCtx->np > 0) {
      LOG_ALLOW(GLOBAL, LOG_INFO, "    - Particle Restart Mode: %s\n", simCtx->particleRestartMode);
    }
    
    // --- Initialize PETSc's internal performance logging stage ---
    ierr = PetscLogDefaultBegin(); CHKERRQ(ierr); // REDUNDANT but safe.
    ierr = PetscMemorySetGetMaximumUsage(); CHKERRQ(ierr);
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finished CreateSimulationContext successfully on rank %d.\n", simCtx->rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PetscMkdirRecursive"
/**
 * @brief Create a directory path recursively using PETSc-compatible error handling.
 */
static PetscErrorCode PetscMkdirRecursive(const char *path)
{
    PetscErrorCode ierr;
    char           tmp_path[PETSC_MAX_PATH_LEN];
    char           *p = NULL;
    size_t         len;
    PetscBool      exists;
 
    PetscFunctionBeginUser;
 
    // Create a mutable copy of the path
    len = strlen(path);
    if (len >= sizeof(tmp_path)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Path is too long to process: %s", path);
    }
    strcpy(tmp_path, path);
 
    // If the path ends with a separator, remove it
    if (tmp_path[len - 1] == '/') {
        tmp_path[len - 1] = 0;
    }
 
    // Iterate through the path, creating each directory level
    for (p = tmp_path + 1; *p; p++) {
        if (*p == '/') {
            *p = 0; // Temporarily terminate the string
 
            // Check if this directory level exists
            ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
            if (!exists) {
                ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
            }
            
            *p = '/'; // Restore the separator
        }
    }
    
    // Create the final, full directory path
    ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
    if (!exists) {
        ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Whether two configured directories denote the same location or nest.
 *
 * @details A bare `strcmp` accepted "output/" against "output", and missed
 *          "output/sub" entirely. Comparison is done on normalized segment sequences so
 *          trailing slashes, "./" prefixes and repeated separators cannot defeat it.
 *
 * @param[in] first  First directory value.
 * @param[in] second Second directory value.
 * @return PETSC_TRUE when one contains the other, or they are equal.
 */
static PetscBool DirectoriesOverlap(const char *first, const char *second)
{
    char        a[PETSC_MAX_PATH_LEN], b[PETSC_MAX_PATH_LEN];
    const char *source[2];
    char       *target[2];
    int         which;
 
    if (!first || !second || !*first || !*second) return PETSC_FALSE;
    source[0] = first;  source[1] = second;
    target[0] = a;      target[1] = b;
 
    for (which = 0; which < 2; ++which) {
        const char *cursor = source[which];
        char       *out    = target[which];
        size_t      used   = 0;
 
        while (*cursor && used + 1 < PETSC_MAX_PATH_LEN) {
            size_t segment;
            while (*cursor == '/') ++cursor;
            segment = strcspn(cursor, "/");
            if (segment == 0) break;
            if (segment == 1 && cursor[0] == '.') { cursor += segment; continue; }
            if (used) out[used++] = '/';
            if (used + segment + 1 >= PETSC_MAX_PATH_LEN) break;
            memcpy(out + used, cursor, segment);
            used += segment;
            cursor += segment;
        }
        out[used] = '\0';
    }
 
    if (!*a || !*b) return PETSC_FALSE;
    if (strcmp(a, b) == 0) return PETSC_TRUE;
    if (strncmp(a, b, strlen(b)) == 0 && a[strlen(b)] == '/') return PETSC_TRUE;
    if (strncmp(b, a, strlen(a)) == 0 && b[strlen(a)] == '/') return PETSC_TRUE;
    return PETSC_FALSE;
}
 
/** @brief Directory names the run tree owns; a log directory must never target one. */
static const char *kReservedRunDirectories[] = {"config", "scheduler", "checkpoints", "visualization", NULL};
 
/**
 * @brief Whether a configured directory name is safe to write to a PETSc options line.
 *
 * @details Whitespace, quotes, and comment markers change how the options file is
 *          tokenized, so a value carrying them would be misread rather than merely
 *          unsafe. This is never waivable: the user cannot consent to ambiguity.
 *
 * @param[in] value Configured directory value.
 * @return PETSC_TRUE when the value is unambiguous.
 */
static PetscBool DirectoryValueIsWellFormed(const char *value)
{
    const char *cursor;
 
    if (!value || value[0] == '\0') return PETSC_FALSE;
    for (cursor = value; *cursor; ++cursor) {
        if (isspace((unsigned char)*cursor)) return PETSC_FALSE;
        if (*cursor == '"' || *cursor == '\'' || *cursor == '#') return PETSC_FALSE;
    }
    return PETSC_TRUE;
}
 
/**
 * @brief Whether a directory's first path segment collides with a reserved run directory.
 *
 * @param[in] value Configured directory value.
 * @return PETSC_TRUE when it collides.
 */
static PetscBool DirectoryHitsReservedName(const char *value)
{
    const char *cursor = value;
 
    if (!value) return PETSC_FALSE;
 
    /* Walk every path segment, skipping separators and "." components. A
       first-segment-only check accepted "./config", ".//config", "output/" and
       "output/sub", each of which deletes a run-owned directory. */
    while (*cursor) {
        size_t segment;
        int    index;
 
        while (*cursor == '/') ++cursor;
        if (!*cursor) break;
        segment = strcspn(cursor, "/");
        if (!(segment == 1 && cursor[0] == '.')) {
            for (index = 0; kReservedRunDirectories[index]; ++index) {
                if (strlen(kReservedRunDirectories[index]) == segment &&
                    strncmp(cursor, kReservedRunDirectories[index], segment) == 0) {
                    return PETSC_TRUE;
                }
            }
        }
        cursor += segment;
    }
    return PETSC_FALSE;
}
 
typedef enum {
    DIR_VERDICT_CONTAINED = 0,      /* Inside the working directory. Always safe.        */
    DIR_VERDICT_MALFORMED,          /* Empty, or unwritable to an options line.          */
    DIR_VERDICT_RESERVED,           /* Collides with a reserved run directory.           */
    DIR_VERDICT_OVERLAP,            /* Nests with the solver output directory.           */
    DIR_VERDICT_RUN_ROOT,           /* Resolves to the working directory itself.         */
    DIR_VERDICT_ANCESTOR,           /* Contains the working directory.                   */
    DIR_VERDICT_RELATIVE_ESCAPE,    /* Relative traversal, or a relative symlink escape. */
    DIR_VERDICT_UNEXPANDED_TILDE,   /* Starts with '~', which nothing expands.           */
    DIR_VERDICT_EXTERNAL_ABSOLUTE,  /* Absolute location outside the tree.               */
    DIR_VERDICT_UNRESOLVABLE        /* Working directory or path could not be resolved.  */
} DirectoryVerdict;
 
/**
 * @brief Lexically normalize a path, resolving "." and ".." textually.
 *
 * @details Used only when the target does not exist yet, so `realpath` cannot resolve
 *          it. Resolving ".." textually is correct here precisely because there is no
 *          file to be a symlink; the existing prefix is resolved separately with
 *          `realpath`, which does follow symlinks.
 *
 * @param[in]  value Path to normalize.
 * @param[out] out   Buffer receiving the normalized path.
 * @param[in]  size  Size of the output buffer.
 * @param[out] stack Caller-provided scratch of at least PETSC_MAX_PATH_LEN bytes.
 * @return PETSC_TRUE when normalization fit in the buffer.
 */
static PetscBool NormalizePathLexically(const char *value, char *out, size_t size,
                                        char *stack)
{
    const char *cursor = value;
    size_t      used = 0;
    PetscBool   absolute = (PetscBool)(value[0] == '/');
 
    stack[0] = '\0';
    while (*cursor) {
        const char *slash;
        size_t      len;
 
        while (*cursor == '/') cursor++;
        if (!*cursor) break;
        slash = strchr(cursor, '/');
        len   = slash ? (size_t)(slash - cursor) : strlen(cursor);
 
        if (len == 1 && cursor[0] == '.') {
            /* Nothing to do. */
        } else if (len == 2 && cursor[0] == '.' && cursor[1] == '.') {
            char *last = strrchr(stack, '/');
            if (last) { *last = '\0'; used = strlen(stack); }
            else if (used) { stack[0] = '\0'; used = 0; }
        } else {
            if (used + len + 2 >= PETSC_MAX_PATH_LEN) return PETSC_FALSE;
            stack[used++] = '/';
            memcpy(stack + used, cursor, len);
            used += len;
            stack[used] = '\0';
        }
        if (!slash) break;
        cursor = slash + 1;
    }
 
    if (absolute) {
        if (strlen(stack) + 1 >= size) return PETSC_FALSE;
        strcpy(out, used ? stack : "/");
    } else {
        const char *relative = used ? stack + 1 : ".";
        if (strlen(relative) + 1 >= size) return PETSC_FALSE;
        strcpy(out, relative);
    }
    return PETSC_TRUE;
}
 
/**
 * @brief Resolve a directory that may not exist yet to an absolute physical path.
 *
 * @details `realpath` fails on a path whose final components have not been created,
 *          which is the normal case for a log directory. The longest existing prefix is
 *          resolved with `realpath`, so symlinked ancestors are followed, and the
 *          remainder is appended lexically.
 *
 * @param[in]  value Configured directory value.
 * @param[in]  cwd     Current working directory, for relative values.
 * @param[out] out     Buffer receiving the resolved absolute path.
 * @param[in]  size    Size of the output buffer.
 * @param[out] scratch Caller-provided scratch of at least 5*PETSC_MAX_PATH_LEN bytes.
 * @return PETSC_TRUE when a resolution was produced.
 */
static PetscBool ResolveDirectoryPhysically(const char *value, const char *cwd,
                                            char *out, size_t size, char *scratch)
{
    /* Five PATH_MAX buffers, on the heap. On the stack they cost 20 KB inside a call
       chain that already carries the caller's three, and that overflowed far enough to
       corrupt the simulation context on rank zero - the guard broke the run it was
       there to protect. A path check must not spend the stack. */
    char *absolute   = scratch;
    char *normalized = scratch + PETSC_MAX_PATH_LEN;
    char *probe      = scratch + 2 * PETSC_MAX_PATH_LEN;
    char *resolved   = scratch + 3 * PETSC_MAX_PATH_LEN;
    char *lexical    = scratch + 4 * PETSC_MAX_PATH_LEN;
 
    if (value[0] == '/') {
        if ((size_t)snprintf(absolute, PETSC_MAX_PATH_LEN, "%s", value) >= PETSC_MAX_PATH_LEN)
            return PETSC_FALSE;
    } else {
        if ((size_t)snprintf(absolute, PETSC_MAX_PATH_LEN, "%s/%s", cwd, value)
            >= PETSC_MAX_PATH_LEN)
            return PETSC_FALSE;
    }
    if (!NormalizePathLexically(absolute, normalized, PETSC_MAX_PATH_LEN, lexical))
        return PETSC_FALSE;
 
    if ((size_t)snprintf(probe, PETSC_MAX_PATH_LEN, "%s", normalized) >= PETSC_MAX_PATH_LEN)
        return PETSC_FALSE;
 
    /* Walk up to the longest existing prefix, then re-append what was trimmed. */
    for (;;) {
        char *last;
        if (realpath(probe, resolved)) {
            size_t consumed = strlen(probe);
            const char *tail = normalized + consumed;
            if ((size_t)snprintf(out, size, "%s%s", resolved, tail) >= size) return PETSC_FALSE;
            /* Re-normalize: `resolved` may be "/" and `tail` may start with "/". */
            return NormalizePathLexically(out, out, size, lexical);
        }
        last = strrchr(probe, '/');
        if (!last) return PETSC_FALSE;
        if (last == probe) { probe[1] = '\0'; }   /* Down to "/" - try once more. */
        else              { *last = '\0'; }
        if (strcmp(probe, "/") == 0 && !realpath(probe, resolved)) return PETSC_FALSE;
    }
}
 
/**
 * @brief Whether `ancestor` is the same directory as `path`, or contains it.
 * @param[in] ancestor Candidate containing directory, absolute and normalized.
 * @param[in] path     Candidate contained directory, absolute and normalized.
 * @return PETSC_TRUE when `ancestor` equals or contains `path`.
 */
static PetscBool PathContainsOrEquals(const char *ancestor, const char *path)
{
    size_t len = strlen(ancestor);
 
    if (strcmp(ancestor, path) == 0) return PETSC_TRUE;
    if (strcmp(ancestor, "/") == 0) return PETSC_TRUE;
    return (PetscBool)(strncmp(path, ancestor, len) == 0 && path[len] == '/');
}
 
/**
 * @brief Classify a configured log directory against the working directory.
 *
 * @details Every non-waivable check - lexical and physical alike - runs before the
 *          caller is allowed to consider an authorization. An earlier version returned
 *          `safe` as soon as it saw an absolute path with authorization set, which let
 *          an authorized run delete its own run directory, an ancestor of it, or the
 *          filesystem root. Classification is now total and the waiver is applied once,
 *          to exactly one verdict.
 *
 * @param[in]  log_dir    Configured log directory.
 * @param[in]  output_dir Configured output directory, for the overlap check.
 * @param[out] reason     Set to a short explanation of the verdict.
 * @return The verdict for this directory.
 */
static DirectoryVerdict ClassifyLogDirectory(const char *log_dir, const char *output_dir,
                                             const char **reason)
{
    /* One heap block for every path buffer this classification needs: three here and
       five for the physical resolution. */
    char            *scratch = NULL;
    char            *cwd, *resolved, *cwd_real;
    DirectoryVerdict verdict;
 
    *reason = "unknown";
 
    /* --- Non-waivable: ambiguous or self-destructive targets. --- */
    if (!log_dir || log_dir[0] == '\0') {
        *reason = "is empty";
        return DIR_VERDICT_MALFORMED;
    }
    if (!DirectoryValueIsWellFormed(log_dir)) {
        *reason = "contains whitespace, a quote, or a comment marker";
        return DIR_VERDICT_MALFORMED;
    }
    if (DirectoryHitsReservedName(log_dir)) {
        *reason = "collides with a reserved run directory";
        return DIR_VERDICT_RESERVED;
    }
    if (DirectoriesOverlap(log_dir, output_dir)) {
        *reason = "overlaps the solver output directory";
        return DIR_VERDICT_OVERLAP;
    }
    /* `~` is refused before anything else looks at the path. Nothing expands it: the
       control file is read by PETSc rather than by a shell, and the code below resolves
       a value not starting with '/' relative to the working directory, so `~/logs`
       names a literal '~' directory inside the run. An earlier version of the launcher's
       physical check expanded it and treated the result as an authorizable external
       location, while nothing that used the value ever did; refusing it here keeps both
       layers describing the same directory. */
    if (log_dir[0] == '~') {
        *reason = "starts with '~', which nothing expands - the options file is read by "
                  "PETSc, not by a shell, so this would name a literal '~' directory "
                  "inside the run. Give a real absolute path. This cannot be overridden";
        return DIR_VERDICT_UNEXPANDED_TILDE;
    }
    /* A relative traversal is refused on the value itself, before any resolution: it
       lands among sibling runs and study members, and no authorization covers it. */
    if (log_dir[0] != '/' &&
        (strcmp(log_dir, "..") == 0 || strncmp(log_dir, "../", 3) == 0 ||
         strstr(log_dir, "/../") != NULL ||
         (strlen(log_dir) >= 3 && strcmp(log_dir + strlen(log_dir) - 3, "/..") == 0))) {
        *reason = "walks above the working directory by relative traversal, which cannot "
                  "be authorized";
        return DIR_VERDICT_RELATIVE_ESCAPE;
    }
 
    if (PetscMalloc1(8 * PETSC_MAX_PATH_LEN, &scratch)) {
        *reason = "cannot be checked because scratch space could not be allocated";
        return DIR_VERDICT_UNRESOLVABLE;
    }
    cwd      = scratch;
    resolved = scratch + PETSC_MAX_PATH_LEN;
    cwd_real = scratch + 2 * PETSC_MAX_PATH_LEN;
 
    verdict = DIR_VERDICT_UNRESOLVABLE;
    if (!getcwd(cwd, PETSC_MAX_PATH_LEN)) {
        *reason = "cannot be checked because the working directory could not be determined";
    } else if (!realpath(cwd, cwd_real)) {
        *reason = "cannot be checked because the working directory could not be resolved";
    } else if (!ResolveDirectoryPhysically(log_dir, cwd_real, resolved, PETSC_MAX_PATH_LEN,
                                           scratch + 3 * PETSC_MAX_PATH_LEN)) {
        *reason = "could not be resolved to a physical path";
    } else if (strcmp(resolved, cwd_real) == 0) {
        /* --- Physical verdicts, decided before authorization is consulted. --- */
        *reason = "resolves to the working directory itself; deleting it would destroy "
                  "the run. This cannot be overridden";
        verdict = DIR_VERDICT_RUN_ROOT;
    } else if (PathContainsOrEquals(resolved, cwd_real)) {
        *reason = "contains the working directory, so deleting it recursively would "
                  "destroy the run and everything beside it. This cannot be overridden";
        verdict = DIR_VERDICT_ANCESTOR;
    } else if (PathContainsOrEquals(cwd_real, resolved)) {
        verdict = DIR_VERDICT_CONTAINED;
    } else if (log_dir[0] == '/') {
        /* Outside the working directory. Only an explicitly absolute value can be
           waived: a relative name that lands outside got there through a symlink,
           which is not a location anybody asked for. */
        *reason = "is an absolute path outside the working directory and no "
                  "-allow_unsafe_log_dir authorization was given";
        verdict = DIR_VERDICT_EXTERNAL_ABSOLUTE;
    } else {
        *reason = "is a relative name that resolves outside the working directory "
                  "through a symlink, which cannot be authorized";
        verdict = DIR_VERDICT_RELATIVE_ESCAPE;
    }
 
    PetscFree(scratch);
    return verdict;
}
 
/**
 * @brief Final safety guard before the runtime deletes its log directory.
 *
 * @details `picurv` validates run directories and re-checks them at submission, but the
 *          solver can be launched directly on a hand-written control file, and a
 *          directory can be replaced with a symlink between validation and launch. This
 *          is the last check before an irreversible recursive delete, so it is
 *          deliberately independent of anything the launcher did, which is what
 *          eliminates the validation-to-execution race.
 *
 *          An explicit `-allow_unsafe_log_dir` authorization waives **exactly one**
 *          verdict: an absolute location outside the working directory that neither is
 *          nor contains it. It never waives the working directory itself, an ancestor
 *          of it, a reserved run directory, an overlap with the output directory, a
 *          relative traversal, a relative symlink escape, or a malformed value. Those
 *          destroy the run itself or corrupt the options file, and no user consent
 *          makes them correct.
 *
 * @param[in]  log_dir    Configured log directory.
 * @param[in]  output_dir Configured output directory, for the overlap check.
 * @param[in]  authorized Whether an explicit unsafe-path authorization was supplied.
 * @param[out] reason     Set to a short explanation when the directory is refused.
 * @return PETSC_TRUE when the directory is safe to remove.
 */
static PetscBool LogDirectoryIsSafeToWipe(const char *log_dir, const char *output_dir,
                                          PetscBool authorized, const char **reason)
{
    DirectoryVerdict verdict = ClassifyLogDirectory(log_dir, output_dir, reason);
 
    if (verdict == DIR_VERDICT_CONTAINED) return PETSC_TRUE;
    if (verdict == DIR_VERDICT_EXTERNAL_ABSOLUTE && authorized) return PETSC_TRUE;
    return PETSC_FALSE;
}
 
#undef __FUNCT__
#define __FUNCT__ "SetupSimulationEnvironment"
/**
 * @brief Internal helper implementation: `SetupSimulationEnvironment()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetupSimulationEnvironment(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscBool      exists;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Setting up simulation environment ---\n");
 
    /* =====================================================================
     *  Phase 1: Check for all required and optional INPUT files.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 1: Verifying input files...\n");
 
    // --- Mandatory Inputs ---
    if (!simCtx->generate_grid) {
        ierr = VerifyPathExistence(simCtx->grid_file, PETSC_FALSE, PETSC_FALSE, "Grid file", &exists); CHKERRQ(ierr);
    }
    for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
        char desc[128];
        ierr = PetscSNPrintf(desc, sizeof(desc), "BCS file #%d", i + 1); CHKERRQ(ierr);
        ierr = VerifyPathExistence(simCtx->bcs_files[i], PETSC_FALSE, PETSC_FALSE, desc, &exists); CHKERRQ(ierr);
    }
 
    // --- Optional Inputs (these produce warnings if missing) ---
    if (simCtx->useCfg) {
        ierr = VerifyPathExistence(simCtx->allowedFile, PETSC_FALSE, PETSC_TRUE, "Whitelist config file", &exists); CHKERRQ(ierr);
    }
    if (simCtx->useProfilingSelectedFuncsCfg) {
        ierr = VerifyPathExistence(simCtx->profilingSelectedFuncsFile, PETSC_FALSE, PETSC_TRUE, "Profiling config file", &exists); CHKERRQ(ierr);
    }
    if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
        ierr = VerifyPathExistence(simCtx->PostprocessingControlFile, PETSC_FALSE, PETSC_TRUE, "Post-processing control file", &exists); CHKERRQ(ierr);
    }
 
 
    /* =====================================================================
     *  Phase 2: Validate directories specific to the execution mode.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 2: Verifying execution mode directories...\n");
    // The data source directory must exist if we intend to load any data from it.
    // This is true if:
    //  1. We are restarting from a previous time step (StartStep > 0), which implies
    //     loading Eulerian fields and/or particle fields.
    //  2. We are starting from t=0 but are explicitly told to load the initial
    //     Eulerian fields from a file (eulerianSource == "load").
    if (simCtx->StartStep > 0 || strcmp(simCtx->eulerianSource,"load")== 0){ // If this is a restart run
        ierr = VerifyPathExistence(simCtx->restart_dir, PETSC_TRUE, PETSC_FALSE, "Restart source directory", &exists); CHKERRQ(ierr);
    }
    if (simCtx->StartStep == 0 && strcmp(simCtx->eulerianSource, "solve") == 0 &&
        simCtx->initialConditionMode == IC_MODE_FILE) {
        ierr = VerifyPathExistence(simCtx->initialConditionDirectory, PETSC_TRUE, PETSC_FALSE,
                                   "Initial-condition source directory", &exists); CHKERRQ(ierr);
    }
    if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
        ierr = VerifyPathExistence(simCtx->pps->source_dir, PETSC_TRUE, PETSC_FALSE, "Post-processing source directory", &exists); CHKERRQ(ierr);
    }
 
    /* =====================================================================
     *  Phase 3: Create and prepare all OUTPUT directories.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 3: Preparing output directories...\n");
 
    if (rank == 0){
      if(simCtx->exec_mode == EXEC_MODE_SOLVER){
        // --- Prepare Log Directory ---
        if (!simCtx->continueMode) {
            // Only wipe logs on fresh runs; continue mode appends to existing logs.
            // Final guard: this delete is recursive and irreversible, so re-check
            // containment here rather than trusting that the launcher validated it.
            {
                PetscBool   unsafe_authorized = PETSC_FALSE;
                const char *refusal = NULL;
                ierr = PetscOptionsGetBool(NULL, NULL, "-allow_unsafe_log_dir",
                                           &unsafe_authorized, NULL); CHKERRQ(ierr);
                if (!LogDirectoryIsSafeToWipe(simCtx->log_dir, simCtx->output_dir,
                                              unsafe_authorized, &refusal)) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
                            "Refusing to delete log directory '%s': it %s. Configure run "
                            "directories through monitor.io.directories.",
                            simCtx->log_dir, refusal);
                }
                if (unsafe_authorized) {
                    LOG_ALLOW(GLOBAL, LOG_WARNING,
                              "Deleting log directory '%s' outside the working directory, "
                              "authorized by -allow_unsafe_log_dir.\n", simCtx->log_dir);
                }
            }
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating/cleaning log directory: %s\n", simCtx->log_dir);
            ierr = PetscRMTree(simCtx->log_dir); // Wipes the directory and its contents
            if (ierr) { /* Ignore file-not-found error, but fail on others */
                PetscError(PETSC_COMM_SELF, __LINE__, __FUNCT__, __FILE__, ierr, PETSC_ERROR_INITIAL, "Could not remove existing log directory '%s'. Check permissions.", simCtx->log_dir);
            }
            ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
        } else {
            // In continue mode, ensure log directory exists but don't wipe it.
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "Continue mode: preserving existing log directory: %s\n", simCtx->log_dir);
            ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
        }
 
        // --- Prepare Output Directory ---
        // The checkpoint coordinator creates its fixed internal hierarchy
        // transactionally when a checkpoint is due.
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying main output directory: %s\n", simCtx->output_dir);
        ierr = PetscTestDirectory(simCtx->output_dir, 'r', &exists); CHKERRQ(ierr);
        if (!exists) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "Output directory not found. Creating: %s\n", simCtx->output_dir);
            ierr = PetscMkdir(simCtx->output_dir); CHKERRQ(ierr);
        }
        // The analysis home is nested (output/analysis/metrics by default), so it
        // needs the recursive create the postprocessor branch already uses; PetscMkdir
        // fails on a path whose parents do not exist yet.
        ierr = PetscMkdirRecursive(simCtx->analysis_dir); CHKERRQ(ierr);
      } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing post-processing output directories ...\n");
 
        PostProcessParams *pps = simCtx->pps;
        char path_buffer[PETSC_MAX_PATH_LEN];
 
        const char *last_slash_euler = strrchr(pps->output_prefix, '/');
        if(last_slash_euler){
          size_t dir_len = last_slash_euler - pps->output_prefix;
          if(dir_len > 0){
            if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing output prefix path is too long.");
            strncpy(path_buffer, pps->output_prefix, dir_len);
            path_buffer[dir_len] = '\0';
 
            ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
            if (!exists){
              LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Eulerian output directory: %s\n", path_buffer);
              ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
            }
          }
        }
 
        // Particle output directory
        if(pps->outputParticles){
          const char  *last_slash_particle = strrchr(pps->particle_output_prefix, '/');
          if(last_slash_particle){
            size_t dir_len = last_slash_particle - pps->particle_output_prefix;
            if(dir_len >  0){
              if(dir_len > sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing particle output prefix path is too long.");
              strncpy(path_buffer, pps->particle_output_prefix, dir_len);
              path_buffer[dir_len] = '\0';
              
              ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
 
              if (!exists){
                LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Particle output directory: %s\n", path_buffer);
                ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
              }
            }
          }
        }
 
        // Statistics output directory
        if(pps->statistics_pipeline[0] != '\0'){
          const char *last_slash_stats = strrchr(pps->statistics_output_prefix, '/');
          if(last_slash_stats){
            size_t dir_len = last_slash_stats - pps->statistics_output_prefix;
            if(dir_len > 0){
              if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing statistics output prefix path is too long.");
              strncpy(path_buffer, pps->statistics_output_prefix, dir_len);
              path_buffer[dir_len] = '\0';
 
              ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
              if (!exists){
                LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Statistics output directory: %s\n", path_buffer);
                ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
              }
            }
          }
        }
      }
    }
    
    // Synchronize all processes before proceeding
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    /* PETSc opens a -info file during PetscInitialize, before the fresh-run wipe above
       removes the log directory it usually lives in. Every later record would then go
       to an unlinked file and the run would end with no info log at all, so each rank
       reopens its file once the directory exists again. Records from initialization
       itself are the only ones lost. */
    if (simCtx->exec_mode == EXEC_MODE_SOLVER && !simCtx->continueMode) {
        char *info_name = NULL;
        FILE *info_file = NULL;
        ierr = PetscInfoGetFile(&info_name, &info_file); CHKERRQ(ierr);
        if (info_name && info_name[0] != '\0') {
            char reopen_name[PETSC_MAX_PATH_LEN];
            ierr = PetscStrncpy(reopen_name, info_name, sizeof(reopen_name)); CHKERRQ(ierr);
            if (info_file && info_file != PETSC_STDOUT) {
                ierr = PetscFClose(PETSC_COMM_SELF, info_file); CHKERRQ(ierr);
            }
            ierr = PetscInfoSetFile(reopen_name, "w"); CHKERRQ(ierr);
        }
        /* PetscInfoGetFile hands back its own copy of the name. */
        ierr = PetscFree(info_name); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Environment setup complete ---\n");
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "AllocateContextHeirarchy"
/**
 * @brief Allocate the user-context objects required by every multigrid level.
 */
static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx;
    PetscInt       nblk = simCtx->block_number;
    PetscBool      found;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Allocating context hierarchy for %d levels and %d blocks...\n", simCtx->mglevels, nblk);
 
    // Store the number of levels in the UserMG struct itself
    usermg->mglevels = simCtx->mglevels;
 
    // --- 1. Allocate the array of MGCtx structs ---
    ierr = PetscMalloc(usermg->mglevels * sizeof(MGCtx), &usermg->mgctx); CHKERRQ(ierr);
    // Zero-initialize to ensure all pointers (especially packer) are NULL
    ierr = PetscMemzero(usermg->mgctx, usermg->mglevels * sizeof(MGCtx)); CHKERRQ(ierr);
    mgctx = usermg->mgctx;
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Allocated MGCtx array of size %d.\n", simCtx->rank, usermg->mglevels);
    
    // --- 2. Parse semi-coarsening options (logic from MG_Initial) ---
    // These flags determine if a dimension is coarsened in the multigrid hierarchy.
    PetscInt *isc, *jsc, *ksc;
    ierr = PetscMalloc3(nblk, &isc, nblk, &jsc, nblk, &ksc); CHKERRQ(ierr);
    // Set defaults to FALSE (full coarsening)
    for (PetscInt i = 0; i < nblk; ++i) {
        isc[i] = 0; jsc[i] = 0; ksc[i] = 0;
    }
 
// Use a temporary variable for the 'count' argument to the parsing function.
    // This protects the original 'nblk' which is needed for the loop bounds.
    PetscInt n_opts_found = nblk;
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_i_semi", isc, &n_opts_found, &found); CHKERRQ(ierr);
 
    n_opts_found = nblk; // Reset the temp variable before the next call
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_j_semi", jsc, &n_opts_found, &found); CHKERRQ(ierr);
 
    n_opts_found = nblk; // Reset the temp variable before the next call
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_k_semi", ksc, &n_opts_found, &found); CHKERRQ(ierr);
 
    // --- 3. Loop over levels and blocks to allocate UserCtx arrays ---
    for (PetscInt level = 0; level < simCtx->mglevels; level++) {
 
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting up MG Level %d...\n", simCtx->rank, level);        
        // Allocate the array of UserCtx structs for this level
        ierr = PetscMalloc(nblk * sizeof(UserCtx), &mgctx[level].user); CHKERRQ(ierr);
        // It's good practice to zero out the memory to avoid uninitialized values
        ierr = PetscMemzero(mgctx[level].user, nblk * sizeof(UserCtx)); CHKERRQ(ierr);
        mgctx[level].thislevel = level;
 
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *currentUser = &mgctx[level].user[bi];
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Initializing UserCtx for Level %d, Block %d.\n", simCtx->rank, level, bi);
        
            // --- CRITICAL STEP: Set the back-pointer to the master context ---
            currentUser->simCtx = simCtx;
 
            // Initialize other per-context values
            currentUser->thislevel = level;
            currentUser->_this = bi; //
            currentUser->mglevels = usermg->mglevels;
 
            // Assign semi-coarsening flags
            currentUser->isc = isc[bi];
            currentUser->jsc = jsc[bi];
            currentUser->ksc = ksc[bi];
 
            // Link to finer/coarser contexts for multigrid operations. Levels are
            // allocated coarsest first, so the finer array does not exist yet when a
            // level is created; each level links itself into the coarser one instead.
            if (level > 0) {
                currentUser->user_c = &mgctx[level-1].user[bi];
                mgctx[level-1].user[bi].user_f = currentUser;
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d:     -> Linked to coarser context (user_c) and back (user_f).\n", simCtx->rank);
            }
        }
    }
 
    // Log a summary of the parsed flags on each rank.
    if (get_log_level() >= LOG_DEBUG && nblk > 0) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Final semi-coarsening configuration view:\n", simCtx->rank);
        for (PetscInt bi = 0; bi < nblk; ++bi) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Block %d: i-semi=%d, j-semi=%d, k-semi=%d\n", simCtx->rank, bi, isc[bi], jsc[bi], ksc[bi]);
        }
    }
 
    // Clean up temporary arrays
    ierr = PetscFree3(isc, jsc, ksc); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Context hierarchy allocation complete.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetupGridAndSolvers"
/**
 * @brief Implementation of \ref SetupGridAndSolvers().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see SetupGridAndSolvers()
 */
PetscErrorCode SetupGridAndSolvers(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Grid and Solvers Setup ---\n");
 
    // Phase 1: Allocate the UserMG and UserCtx hierarchy
    ierr = AllocateContextHierarchy(simCtx); CHKERRQ(ierr);
 
    ierr = DefineAllGridDimensions(simCtx); CHKERRQ(ierr);
    ierr = InitializeAllGridDMs(simCtx); CHKERRQ(ierr);
    ierr = AssignAllGridCoordinates(simCtx); CHKERRQ(ierr);
    ierr = CreateAndInitializeAllVectors(simCtx); CHKERRQ(ierr);
 
    // NOTE: CalculateAllGridMetrics is now called inside SetupBoundaryConditions (not here) to ensure:
    // 1. Boundary condition configuration data (boundary_faces) is available for periodic BC corrections
    // 2. Computed metrics are available for inlet/outlet area calculations
    // This resolves the circular dependency between BC setup and metric calculations.
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Grid and Solvers Setup Complete ---\n");
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "CreateAndInitializeAllVectors"
/**
 * @brief Internal helper implementation: `CreateAndInitializeAllVectors()`.
 * @details Local to this translation unit.
 */
PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx = usermg->mgctx;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Creating and initializing all simulation vectors...\n");
        
    for (PetscInt level  = usermg->mglevels-1; level >=0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *user = &mgctx[level].user[bi];
 
            if(!user->da || !user->fda) {
              SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMs not properly initialized in UserCtx before vector creation.");
            }
 
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating vectors for level %d, block %d\n", simCtx->rank, level, bi);
 
        // --- Group A: Primary Flow Fields (Global and Local) ---
            // These are the core solution variables.
            ierr = DMCreateGlobalVector(user->fda, &user->Ucont); CHKERRQ(ierr); ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->Ucat);  CHKERRQ(ierr); ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->P);     CHKERRQ(ierr); ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->Nvert); CHKERRQ(ierr); ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
 
            ierr = DMCreateLocalVector(user->fda, &user->lUcont); CHKERRQ(ierr); ierr = VecSet(user->lUcont, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda, &user->lUcat);  CHKERRQ(ierr); ierr = VecSet(user->lUcat, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lP);     CHKERRQ(ierr); ierr = VecSet(user->lP, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lNvert); CHKERRQ(ierr); ierr = VecSet(user->lNvert, 0.0); CHKERRQ(ierr);
 
            // -- Group A2: Derived Flow Fields (Global and Local) ---
            ierr = VecDuplicate(user->P,&user->Diffusivity); CHKERRQ(ierr); ierr = VecSet(user->Diffusivity, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lP,&user->lDiffusivity); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivity, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Ucat,&user->DiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->DiffusivityGradient, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lUcat,&user->lDiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivityGradient, 0.0); CHKERRQ(ierr);
 
            // -- Group B: Solver Work Vectors (Global and Local) ---
            ierr = VecDuplicate(user->P, &user->Phi);       CHKERRQ(ierr); ierr = VecSet(user->Phi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lP, &user->lPhi);       CHKERRQ(ierr); ierr = VecSet(user->lPhi, 0.0); CHKERRQ(ierr);
 
            // --- Group C: Time-Stepping & Workspace Fields (Finest Level Only) ---
            if (level == usermg->mglevels - 1) {
                ierr = VecDuplicate(user->Ucont, &user->Ucont_o);   CHKERRQ(ierr); ierr = VecSet(user->Ucont_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->Ucont, &user->Ucont_rm1); CHKERRQ(ierr); ierr = VecSet(user->Ucont_rm1, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->Ucat,  &user->Ucat_o);    CHKERRQ(ierr); ierr = VecSet(user->Ucat_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->P,     &user->P_o);       CHKERRQ(ierr); ierr = VecSet(user->P_o, 0.0); CHKERRQ(ierr);
                  ierr = VecDuplicate(user->lUcont, &user->lUcont_o);   CHKERRQ(ierr); ierr = VecSet(user->lUcont_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->lUcont, &user->lUcont_rm1); CHKERRQ(ierr); ierr = VecSet(user->lUcont_rm1, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNvert_o); CHKERRQ(ierr); ierr = VecSet(user->lNvert_o, 0.0); CHKERRQ(ierr);
                    ierr = VecDuplicate(user->Nvert, &user->Nvert_o); CHKERRQ(ierr); ierr = VecSet(user->Nvert_o, 0.0); CHKERRQ(ierr);
              }
 
        // --- Group D: Grid Metrics (Face-Centered) ---
            ierr = DMCreateGlobalVector(user->fda, &user->Csi); CHKERRQ(ierr); ierr = VecSet(user->Csi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->Eta);         CHKERRQ(ierr); ierr = VecSet(user->Eta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->Zet);         CHKERRQ(ierr); ierr = VecSet(user->Zet, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->Aj);  CHKERRQ(ierr); ierr = VecSet(user->Aj, 0.0); CHKERRQ(ierr);
            
            ierr = DMCreateLocalVector(user->fda, &user->lCsi); CHKERRQ(ierr); ierr = VecSet(user->lCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lEta);       CHKERRQ(ierr); ierr = VecSet(user->lEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lZet);       CHKERRQ(ierr); ierr = VecSet(user->lZet, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lAj);  CHKERRQ(ierr); ierr = VecSet(user->lAj, 0.0); CHKERRQ(ierr);
         
        
       // --- Group E: Grid Metrics (Face-Centered) ---
            // Vector metrics are duplicated from Csi (DOF=3, fda-based)
            ierr = VecDuplicate(user->Csi, &user->ICsi); CHKERRQ(ierr); ierr = VecSet(user->ICsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->IEta); CHKERRQ(ierr); ierr = VecSet(user->IEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->IZet); CHKERRQ(ierr); ierr = VecSet(user->IZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JCsi); CHKERRQ(ierr); ierr = VecSet(user->JCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JEta); CHKERRQ(ierr); ierr = VecSet(user->JEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JZet); CHKERRQ(ierr); ierr = VecSet(user->JZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KCsi); CHKERRQ(ierr); ierr = VecSet(user->KCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KEta); CHKERRQ(ierr); ierr = VecSet(user->KEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KZet); CHKERRQ(ierr); ierr = VecSet(user->KZet, 0.0); CHKERRQ(ierr);
            // Scalar metrics are duplicated from Aj (DOF=1, da-based)
            ierr = VecDuplicate(user->Aj, &user->IAj); CHKERRQ(ierr); ierr = VecSet(user->IAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Aj, &user->JAj); CHKERRQ(ierr); ierr = VecSet(user->JAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Aj, &user->KAj); CHKERRQ(ierr); ierr = VecSet(user->KAj, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->lCsi, &user->lICsi); CHKERRQ(ierr); ierr = VecSet(user->lICsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lIEta); CHKERRQ(ierr); ierr = VecSet(user->lIEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lIZet); CHKERRQ(ierr); ierr = VecSet(user->lIZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJCsi); CHKERRQ(ierr); ierr = VecSet(user->lJCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJEta); CHKERRQ(ierr); ierr = VecSet(user->lJEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJZet); CHKERRQ(ierr); ierr = VecSet(user->lJZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKCsi); CHKERRQ(ierr); ierr = VecSet(user->lKCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKEta); CHKERRQ(ierr); ierr = VecSet(user->lKEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKZet); CHKERRQ(ierr); ierr = VecSet(user->lKZet, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->lAj, &user->lIAj); CHKERRQ(ierr); ierr = VecSet(user->lIAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lAj, &user->lJAj); CHKERRQ(ierr); ierr = VecSet(user->lJAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lAj, &user->lKAj); CHKERRQ(ierr); ierr = VecSet(user->lKAj, 0.0); CHKERRQ(ierr);
        
        // --- Group F: Cell/Face Center Coordinates and Grid Spacing ---
            ierr = DMCreateGlobalVector(user->fda, &user->Cent); CHKERRQ(ierr); ierr = VecSet(user->Cent, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda, &user->lCent); CHKERRQ(ierr); ierr = VecSet(user->lCent, 0.0); CHKERRQ(ierr);
            
            ierr = VecDuplicate(user->Cent, &user->GridSpace); CHKERRQ(ierr); ierr = VecSet(user->GridSpace, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lGridSpace); CHKERRQ(ierr); ierr = VecSet(user->lGridSpace, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->Cent, &user->Centx); CHKERRQ(ierr); ierr = VecSet(user->Centx, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Cent, &user->Centy); CHKERRQ(ierr); ierr = VecSet(user->Centy, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Cent, &user->Centz); CHKERRQ(ierr); ierr = VecSet(user->Centz, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCentx); CHKERRQ(ierr); ierr = VecSet(user->lCentx, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCenty); CHKERRQ(ierr); ierr = VecSet(user->lCenty, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCentz); CHKERRQ(ierr); ierr = VecSet(user->lCentz, 0.0); CHKERRQ(ierr);
 
        if(level == usermg->mglevels -1){
        // --- Group G: Turbulence Models (Finest Level Only) ---
            if (simCtx->les) {
                ierr = DMCreateGlobalVector(user->da, &user->Nu_t); CHKERRQ(ierr); ierr = VecSet(user->Nu_t, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNu_t); CHKERRQ(ierr); ierr = VecSet(user->lNu_t, 0.0); CHKERRQ(ierr);
                LOG_ALLOW(GLOBAL, LOG_DEBUG, "Turbulence viscosity (Nu_t) vectors created for the LES model.\n");
                // Only the dynamic model carries a coefficient field. The constant model
                // reads its coefficient straight from configuration, so allocating,
                // synchronizing, and checkpointing a field of one repeated number would
                // buy nothing.
                if(simCtx->les == DYNAMIC_SMAGORINSKY){
                ierr = DMCreateGlobalVector(user->da,&user->CS); CHKERRQ(ierr); ierr = VecSet(user->CS,0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da,&user->lCs); CHKERRQ(ierr); ierr = VecSet(user->lCs,0.0); CHKERRQ(ierr);
                LOG_ALLOW(GLOBAL, LOG_DEBUG, "Dynamic Smagorinsky coefficient (CS) vectors created.\n");
                }
 
                  // Add K_Omega etc. here as needed
 
                // Note: Add any other vectors from the legacy MG_Initial here as needed.
                // For example: Rhs, Forcing, turbulence Vecs (K_Omega, Nu_t)...
        
            }
 
            /* A wall model is configured independently of LES and RANS - `wall_function`
               is their sibling in the schema, not their child - so its storage cannot sit
               inside their block. It did, which left `ApplyWallFunction` dereferencing a
               null vector on any case that enabled a wall model without one of them. */
            if (simCtx->wallfunction) {
                /* Scalar per cell, so it lives on `da`. It was previously created from
                   `fda`, whose dof is 3, while every reader opens it through `da` - a
                   pairing PETSc rejects outright. */
                ierr = DMCreateGlobalVector(user->da, &user->Friction_Velocity); CHKERRQ(ierr);
                ierr = VecSet(user->Friction_Velocity, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lFriction_Velocity); CHKERRQ(ierr);
                ierr = VecSet(user->lFriction_Velocity, 0.0); CHKERRQ(ierr);
                ierr = DMCreateGlobalVector(user->da, &user->Nu_Wall); CHKERRQ(ierr);
                ierr = VecSet(user->Nu_Wall, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNu_Wall); CHKERRQ(ierr);
                ierr = VecSet(user->lNu_Wall, 0.0); CHKERRQ(ierr);
            }
        // --- Group H: Particle Methods    
        if(simCtx->np>0){
          ierr = DMCreateGlobalVector(user->da,&user->ParticleCount); CHKERRQ(ierr); ierr = VecSet(user->ParticleCount,0.0); CHKERRQ(ierr);
        ierr = DMCreateLocalVector(user->da,&user->lParticleCount); CHKERRQ(ierr); ierr = VecSet(user->lParticleCount,0.0); CHKERRQ(ierr);
        // Scalar field to hold particle scalar property (e.g., temperature, concentration)
          ierr = DMCreateGlobalVector(user->da,&user->Psi); CHKERRQ(ierr); ierr = VecSet(user->Psi,0.0); CHKERRQ(ierr);
        ierr = DMCreateLocalVector(user->da,&user->lPsi); CHKERRQ(ierr); ierr = VecSet(user->lPsi,0.0); CHKERRQ(ierr);
          LOG_ALLOW(GLOBAL,LOG_DEBUG,"ParticleCount & Scalar(Psi) created for %d particles.\n",simCtx->np);
          }
        }
        // --- Group I: Boundary Condition vectors ---
        ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Ubcs); CHKERRQ(ierr);
        ierr = VecSet(user->Bcs.Ubcs, 0.0); CHKERRQ(ierr);
        ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Uch); CHKERRQ(ierr);
        ierr = VecSet(user->Bcs.Uch, 0.0); CHKERRQ(ierr);
        
        // --- Group L: Field-Statistics Accumulators (Finest Level Only) ---
        // Config-counted, like the convergence-state reference fields: the window
        // count is resolved before this factory runs, and each accumulator is
        // duplicated from a vector created above so it inherits DM and layout.
        if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
            ierr = PetscCalloc1((size_t)simCtx->fieldStatisticsWindowCount,
                                &user->fieldStatisticsStorage); CHKERRQ(ierr);
            for (PetscInt w = 0; w < simCtx->fieldStatisticsWindowCount; ++w) {
                ierr = PicurvWindowStorageCreate(user,
                           &simCtx->fieldStatisticsWindows[w].definition,
                           &user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
            }
            LOG_ALLOW(LOCAL, LOG_INFO,
                      "Allocated accumulators for %d statistics window(s).\n",
                      simCtx->fieldStatisticsWindowCount);
        }
 
        // --- Group K: Corner-Staging Workspace (Finest Level Only) ---
        // One pair per block size. Previously these were created lazily inside
        // the interpolation routine and rebuilt whenever the block size changed;
        // allocating both here removes that churn and the cached-size check.
        if (level == usermg->mglevels - 1) {
            ierr = DMCreateGlobalVector(user->da,  &user->CellScalarAtCorner);  CHKERRQ(ierr);
            ierr = VecSet(user->CellScalarAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,   &user->lCellScalarAtCorner); CHKERRQ(ierr);
            ierr = VecSet(user->lCellScalarAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->CellVectorAtCorner);  CHKERRQ(ierr);
            ierr = VecSet(user->CellVectorAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda,  &user->lCellVectorAtCorner); CHKERRQ(ierr);
            ierr = VecSet(user->lCellVectorAtCorner, 0.0); CHKERRQ(ierr);
        }
 
        /* --- Group N: Statistics Output Staging (Finest Level Only) ---
         * The same kind of object as the corner workspace: a named buffer the
         * shared kernels can address, rather than simulation state. It exists only
         * to carry derived statistics to the nodal and VTK paths, so a run without
         * statistics allocates none of it. */
        if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
            ierr = DMCreateGlobalVector(user->da,  &user->PostScalar);       CHKERRQ(ierr);
            ierr = VecSet(user->PostScalar, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,   &user->lPostScalar);      CHKERRQ(ierr);
            ierr = VecSet(user->lPostScalar, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->PostScalarNodal);  CHKERRQ(ierr);
            ierr = VecSet(user->PostScalarNodal, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->PostVector);       CHKERRQ(ierr);
            ierr = VecSet(user->PostVector, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda,  &user->lPostVector);      CHKERRQ(ierr);
            ierr = VecSet(user->lPostVector, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->PostVectorNodal);  CHKERRQ(ierr);
            ierr = VecSet(user->PostVectorNodal, 0.0); CHKERRQ(ierr);
        }
 
      if(level == usermg->mglevels - 1){
        if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Post-processor mode detected. Allocating derived field vectors.\n");
 
                ierr = VecDuplicate(user->P, &user->P_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->P_nodal, 0.0); CHKERRQ(ierr);
 
                ierr = VecDuplicate(user->Ucat, &user->Ucat_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Ucat_nodal, 0.0); CHKERRQ(ierr);
                
                ierr = VecDuplicate(user->P, &user->Qcrit); CHKERRQ(ierr);
                ierr = VecSet(user->Qcrit, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lQcrit); CHKERRQ(ierr);
                ierr = VecSet(user->lQcrit, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->P, &user->Qcrit_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Qcrit_nodal, 0.0); CHKERRQ(ierr);
 
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vectors P_nodal, Ucat_nodal, and Qcrit created.\n");
 
                if(simCtx->np>0){
                ierr = VecDuplicate(user->Psi, &user->Psi_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Psi_nodal, 0.0); CHKERRQ(ierr);
 
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vector Psi_nodal created for particle scalar property.\n");
 
                }
        }else{
                user->P_nodal = NULL;
                user->Ucat_nodal = NULL;
                user->Qcrit = NULL;
                user->lQcrit = NULL;
                user->Qcrit_nodal = NULL;
                user->Psi_nodal = NULL;
        }
      }
    
  }
}
 
    /* Config-counted vectors belong to this factory too. The convergence state
     * allocates a Vec array whose length comes from configuration, so it is
     * invoked here rather than by the caller, keeping every vector this run
     * owns created in one place. It depends only on option-parsed configuration
     * and on the vectors created above. */
    ierr = InitializeSolutionConvergenceState(simCtx); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "RepairPeriodicNormalFaceGhosts"
/**
 * @brief Repairs the adjacent normal ghost layer for periodic face-staggered data.
 *
 * PETSc wraps every component with cell-style indexing. A face family instead
 * needs its adjacent normal ghost shifted by one additional physical face. With
 * the width-three periodic DMDA, the required value is available in the deeper
 * PETSc ghost at -3 or n+2. Tangential ghosts retain PETSc's native wraparound.
 */
static PetscErrorCode RepairPeriodicNormalFaceGhosts(UserCtx *user, DM dm, Vec local_vec,
                                                      PetscInt dof, char face_direction,
                                                      PetscBool component_staggered)
{
    DMDALocalInfo info;
    PetscInt xs, xe, ys, ye, zs, ze;
    PetscInt gxs, gxe, gys, gye, gzs, gze;
    PetscInt mx, my, mz;
 
    PetscFunctionBeginUser;
    if (!face_direction && !component_staggered) PetscFunctionReturn(0);
 
    PetscCall(DMDAGetLocalInfo(dm, &info));
    xs = info.xs; xe = info.xs + info.xm;
    ys = info.ys; ye = info.ys + info.ym;
    zs = info.zs; ze = info.zs + info.zm;
    gxs = info.gxs; gxe = info.gxs + info.gxm;
    gys = info.gys; gye = info.gys + info.gym;
    gzs = info.gzs; gze = info.gzs + info.gzm;
    mx = info.mx; my = info.my; mz = info.mz;
 
    if (component_staggered) {
        Cmpnts ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
            PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                array[k][j][-1].x = array[k][j][-3].x;
        }
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
            PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                array[k][j][mx].x = array[k][j][mx + 2].x;
        }
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
            PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                array[k][-1][i].y = array[k][-3][i].y;
        }
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
            PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                array[k][my][i].y = array[k][my + 2][i].y;
        }
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
            PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                array[-1][j][i].z = array[-3][j][i].z;
        }
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
            PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                array[mz][j][i].z = array[mz + 2][j][i].z;
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
        PetscFunctionReturn(0);
    }
 
    if (dof == 1) {
        PetscReal ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (face_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
                PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][-1] = array[k][j][-3];
            }
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
                PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][mx] = array[k][j][mx + 2];
            }
        } else if (face_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
                PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][-1][i] = array[k][-3][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
                PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][my][i] = array[k][my + 2][i];
            }
        } else if (face_direction == 'k') {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
                PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[-1][j][i] = array[-3][j][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
                PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[mz][j][i] = array[mz + 2][j][i];
            }
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
    } else {
        Cmpnts ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (face_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
                PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][-1] = array[k][j][-3];
            }
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
                PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][mx] = array[k][j][mx + 2];
            }
        } else if (face_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
                PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][-1][i] = array[k][-3][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
                PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][my][i] = array[k][my + 2][i];
            }
        } else if (face_direction == 'k') {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
                PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[-1][j][i] = array[-3][j][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
                PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[mz][j][i] = array[mz + 2][j][i];
            }
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "UpdateLocalGhosts"
/**
 * @brief Updates a catalogued field's local ghost representation.
 * @details Resolves the typed field view, performs the existing PETSc scatter,
 *          and applies the established periodic normal-face repair when needed.
 */
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    FieldView      field_view;
    const char    *field_name;
    Vec            globalVec;
    Vec            localVec;
    DM             dm;
    PetscInt       dof;
    char           face_direction = '\0';
    PetscBool      component_staggered = PETSC_FALSE;
 
    PetscFunctionBeginUser; // Use User version for application code
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = FieldGetView(user, field_id, &field_view); CHKERRQ(ierr);
    PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_GHOST_UPDATE) != 0u,
               PETSC_COMM_SELF, PETSC_ERR_SUP,
               "Field '%s' does not support ghost updates.",
               field_view.descriptor->canonical_name);
    field_name = field_view.descriptor->canonical_name;
    globalVec = field_view.global_vec;
    localVec = field_view.local_vec;
    dm = field_view.dm;
    dof = field_view.descriptor->dof;
 
    switch (field_view.descriptor->sync_class) {
        case FIELD_SYNC_I_FACE:
            face_direction = 'i';
            break;
        case FIELD_SYNC_J_FACE:
            face_direction = 'j';
            break;
        case FIELD_SYNC_K_FACE:
            face_direction = 'k';
            break;
        case FIELD_SYNC_COMPONENT_STAGGERED:
            component_staggered = PETSC_TRUE;
            break;
        case FIELD_SYNC_STANDARD:
            break;
        default:
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
                    "Field '%s' has an invalid ghost synchronization class.", field_name);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, field_name);
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
              rank, field_name, (void*)dm, (void*)globalVec, (void*)localVec);
 
    // --- 3. Optional Debugging: Norm Before Update ---
    // Use your logging convention check
    // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
    if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){
        PetscReal norm_global_before;
        ierr = VecNorm(globalVec, NORM_INFINITY, &norm_global_before); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Global) BEFORE Ghost Update: %g\n", field_name, norm_global_before);
        // Optional: Norm of local vector before update (might contain old ghost values)
        // PetscReal norm_local_before;
        // ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_before); CHKERRQ(ierr);
        // LOG_ALLOW(GLOBAL, LOG_DEBUG,"Max norm '%s' (Local) BEFORE Ghost Update: %g\n", field_name, norm_local_before);
    }
 
    // --- 4. Perform the Global-to-Local Transfer (Ghost Update) ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calling DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
    ierr = DMGlobalToLocalBegin(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
    ierr = RepairPeriodicNormalFaceGhosts(user, dm, localVec, dof, face_direction,
                                          component_staggered); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Completed DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
 
    // --- 5. Optional Debugging: Norm After Update ---
    // Use your logging convention check
    // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
    if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){ // Using your specific check
        PetscReal norm_local_after;
        ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_after); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Local) AFTER Ghost Update: %g\n", field_name, norm_local_after);
 
        // --- 6. Optional Debugging: Specific Point Checks (Example for Ucat on Rank 0/1) ---
        //    (Keep this conditional if it's only for specific debug scenarios)
        if (field_id == FIELD_ID_UCAT) { // Only do detailed checks for Ucat for now
           PetscMPIInt rank_test;
           MPI_Comm_rank(PETSC_COMM_WORLD, &rank_test);
 
           // Get Local Info needed for indexing checks
           DMDALocalInfo info_check;
           ierr = DMDAGetLocalInfo(dm, &info_check); CHKERRQ(ierr); // Use the correct dm
 
           // Buffer for array pointer
           Cmpnts ***lUcat_arr_test = NULL;
           PetscErrorCode ierr_test = 0;
 
           LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Testing '%s' access immediately after ghost update...\n", rank_test, field_name);
           ierr_test = DMDAVecGetArrayDOFRead(dm, localVec, &lUcat_arr_test); // Use correct dm and localVec
 
           if (ierr_test) {
               LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d getting '%s' array after ghost update!\n", rank_test, ierr_test, field_name);
           } else if (!lUcat_arr_test) {
                LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, field_name);
           }
           else {
               // Check owned interior point (e.g., first interior point)
               PetscInt k_int = info_check.zs + (info_check.zm > 1 ? 1 : 0); // Global k index (at least zs+1 if possible)
               PetscInt j_int = info_check.ys + (info_check.ym > 1 ? 1 : 0); // Global j index
               PetscInt i_int = info_check.xs + (info_check.xm > 1 ? 1 : 0); // Global i index
                // Ensure indices are within global bounds if domain is very small
               //if (k_int >= info_check.mz-1) k_int = info_check.mz-2; if (k_int < 1) k_int = 1;
               //if (j_int >= info_check.my-1) j_int = info_check.my-2; if (j_int < 1) j_int = 1;
           // if (i_int >= info_check.mx-1) i_int = info_check.mx-2; if (i_int < 1) i_int = 1;
           // clamp k_int to [1 .. mz-2] 
           if (k_int >= info_check.mz - 1) {
         k_int = info_check.mz - 2;
           }
           if (k_int < 1) {
         k_int = 1;
           }
 
           // clamp j_int to [1 .. my-2] 
           if (j_int >= info_check.my - 1) {
         j_int = info_check.my - 2;
           }
           if (j_int < 1) {
         j_int = 1;
           }
 
           // clamp i_int to [1 .. mx-2]
           if (i_int >= info_check.mx - 1) {
         i_int = info_check.mx - 2;
           }
           if (i_int < 1) {
         i_int = 1;
           }
 
               // Only attempt read if indices are actually owned (relevant for multi-rank)
               if (k_int >= info_check.zs && k_int < info_check.zs + info_check.zm &&
                   j_int >= info_check.ys && j_int < info_check.ys + info_check.ym &&
                   i_int >= info_check.xs && i_int < info_check.xs + info_check.xm)
               {
                  LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read OWNED INTERIOR [%d][%d][%d] (Global)\n", rank_test, k_int, j_int, i_int);
                  Cmpnts test_val_owned_interior = lUcat_arr_test[k_int][j_int][i_int];
                  LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned interior: x=%g\n", rank_test, test_val_owned_interior.x);
               } else {
                  LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping interior test read for non-owned index [%d][%d][%d].\n", rank_test, k_int, j_int, i_int);
               }
 
 
               // Check owned boundary point (e.g., first owned point)
               PetscInt k_bnd = info_check.zs; // Global k index
               PetscInt j_bnd = info_check.ys; // Global j index
               PetscInt i_bnd = info_check.xs; // Global i index
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read OWNED BOUNDARY [%d][%d][%d] (Global)\n", rank_test, k_bnd, j_bnd, i_bnd);
               Cmpnts test_val_owned_boundary = lUcat_arr_test[k_bnd][j_bnd][i_bnd];
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned boundary: x=%g\n", rank_test, test_val_owned_boundary.x);
 
 
               // Check ghost point (e.g., one layer below in k, if applicable)
               if (info_check.zs > 0) { // Only if there's a rank below
                   PetscInt k_ghost = info_check.zs - 1;
                   PetscInt j_ghost = info_check.ys; // Use start of owned y, simple example
                   PetscInt i_ghost = info_check.xs; // Use start of owned x, simple example
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read GHOST [%d][%d][%d] (Global)\n", rank_test, k_ghost, j_ghost, i_ghost);
                   Cmpnts test_val_ghost = lUcat_arr_test[k_ghost][j_ghost][i_ghost];
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading ghost: x=%g\n", rank_test, test_val_ghost.x);
               } else {
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping ghost test read (zs=0).\n", rank_test);
               }
 
               // Restore the array
               ierr_test = DMDAVecRestoreArrayDOFRead(dm, localVec, &lUcat_arr_test);
               if(ierr_test){ LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d restoring '%s' array after test read!\n", rank_test, ierr_test, field_name); }
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, field_name);
           }
        } // end if Ucat
    } // end debug logging check
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, field_name);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetupBoundaryConditions"
/**
 * @brief Internal helper implementation: `SetupBoundaryConditions()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetupBoundaryConditions(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "--- Setting up Boundary Conditions ---\n");
    // --- Phase 1: Parse and initialize BC configuration for all blocks ---
    LOG_ALLOW(GLOBAL,LOG_INFO,"Parsing BC configuration files and initializing boundary condition data structures.\n");
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        LOG_ALLOW(GLOBAL,LOG_DEBUG, "  -> Processing Block %d:\n", bi);
 
    // --- Generate the filename for the current block ---
    const char *current_bc_filename = simCtx->bcs_files[bi];
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"  -> Processing Block %d using config file '%s'\n", bi, current_bc_filename);
       // This will populate user_finest[bi].boundary_faces
 
        //ierr = ParseAllBoundaryConditions(&user_finest[bi],current_bc_filename); CHKERRQ(ierr);
 
        ierr = BoundarySystem_Initialize(&user_finest[bi], current_bc_filename); CHKERRQ(ierr);
    }
 
    // Propogate BC Configuration to coarser levels.
    ierr = PropagateBoundaryConfigToCoarserLevels(simCtx); CHKERRQ(ierr);
 
    // Validate the geometric contract before any metric consumes periodic geometry.
    for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
        UserCtx *level_users = simCtx->usermg.mgctx[level].user;
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            ierr = ValidatePeriodicGeometry(&level_users[bi]); CHKERRQ(ierr);
        }
    }
 
    // --- Calculate Grid Metrics (requires BC configuration) ---
    // NOTE: This MUST be called here (after BC initialization but before inlet/outlet calculations) because:
    // 1. Periodic BC corrections in metric calculations need boundary_faces data to be populated
    // 2. Inlet/Outlet area calculations (below) require computed metrics (Csi, Eta, Zet) to be available
    // Previously this was in SetupGridAndSolvers, but that caused metrics to be computed without BC info.
    LOG_ALLOW(GLOBAL,LOG_INFO,"Computing grid metrics with boundary condition information.\n");
    ierr = CalculateAllGridMetrics(simCtx); CHKERRQ(ierr);
 
    // --- Phase 2: Calculate inlet/outlet properties (requires computed metrics) ---
    LOG_ALLOW(GLOBAL,LOG_INFO,"Calculating inlet and outlet face properties.\n");
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        // Call the function to calculate the center of the inlet face & the inlet area, which may be used to calculate Boundary values.
        ierr = CalculateInletProperties(&user_finest[bi]); CHKERRQ(ierr);
 
        // Call the function to calculate the center of the outlet face & the outlet area, which may be used to calculate Boundary values.
        ierr = CalculateOutletProperties(&user_finest[bi]); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "--- Boundary Conditions setup complete ---\n");   
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Internal helper implementation: `Allocate3DArrayScalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode Allocate3DArrayScalar(PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
{
  PetscErrorCode ierr;
  PetscReal      ***data;
  PetscReal      *dataContiguous;
  PetscInt       k, j;
 
  PetscFunctionBegin;
  /* Step 1: Allocate memory for an array of nz layer pointers (zero-initialized) */
  ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
 
  /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
  ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
  for (k = 1; k < nz; k++) {
    data[k] = data[0] + k * ny;
  }
 
  /* Step 3: Allocate one contiguous block for all data elements (nz*ny*nx) */
  ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
 
  /* Build the 3D pointer structure: each row pointer gets the correct segment of data */
  for (k = 0; k < nz; k++) {
    for (j = 0; j < ny; j++) {
      data[k][j] = dataContiguous + (k * ny + j) * nx;
      /* Memory is already zeroed by PetscCalloc1, so no manual initialization is needed */
    }
  }
  *array = data;
  PetscFunctionReturn(0);
}
 
/**
 * @brief Internal helper implementation: `Deallocate3DArrayScalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode Deallocate3DArrayScalar(PetscReal ***array, PetscInt nz, PetscInt ny)
{
  PetscErrorCode ierr;
  (void)nz;
  (void)ny;
 
  PetscFunctionBegin;
   if (!array || !array[0] || !array[0][0] ) { // Added more robust check
      LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayScalar called with potentially unallocated or NULL array.\n");
       if (array) {
           if (array[0]) { // Check if row pointers might exist
               // Cannot safely access array[0][0] if array[0] might be invalid/freed
               // Standard deallocation below assumes valid pointers.
                ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers if they exist
           }
            ierr = PetscFree(array); CHKERRQ(ierr); // Free layer pointers if they exist
       }
       PetscFunctionReturn(0);
   }
 
  // --- Standard Deallocation (assuming valid allocation) ---
 
  /* 1. Free the contiguous block of PetscReal values.
     The starting address was stored in array[0][0]. */
  ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
 
  /* 2. Free the contiguous block of row pointers.
     The starting address was stored in array[0]. */
  ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
 
  /* 3. Free the layer pointer array.
     The starting address is 'array' itself. */
  ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
 
  PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref Allocate3DArrayVector().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see Allocate3DArrayVector()
 */
PetscErrorCode Allocate3DArrayVector(Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
{
  PetscErrorCode ierr;
  Cmpnts         ***data;
  Cmpnts         *dataContiguous;
  PetscInt       k, j;
  PetscMPIInt    rank;
 
  PetscFunctionBegin;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  
  /* Step 1: Allocate memory for nz layer pointers (zeroed) */
  ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
 
  LOG_ALLOW(LOCAL,LOG_DEBUG," [Rank %d] memory allocated for outermost layer (%d k-layer pointers).\n",rank,nz);
  
  /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
  ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
  for (k = 1; k < nz; k++) {
    data[k] = data[0] + k * ny;
  }
 
  LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] memory allocated for %dx%d row pointers.\n",rank,nz,ny);
  
  /* Step 3: Allocate one contiguous block for nz*ny*nx Cmpnts structures (zeroed) */
  ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] memory allocated for contigous block of %dx%dx%d Cmpnts structures).\n",rank,nz,ny,nx); 
 
  /* Build the 3D pointer structure for vector data */
  for (k = 0; k < nz; k++) {
    for (j = 0; j < ny; j++) {
      data[k][j] = dataContiguous + (k * ny + j) * nx;
      /* The PetscCalloc1 call has already initialized each Cmpnts to zero. */
    }
  }
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] 3D pointer structure for vector data created. \n",rank);
  
  *array = data;
  PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref Deallocate3DArrayVector().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see Deallocate3DArrayVector()
 */
 PetscErrorCode Deallocate3DArrayVector(Cmpnts ***array, PetscInt nz, PetscInt ny)
{
  PetscErrorCode ierr;
  (void)nz;
  (void)ny;
 
  PetscFunctionBegin;
  // If array is NULL or hasn't been allocated properly, just return.
  if (!array || !array[0] || !array[0][0] ) {
      LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayVector called with potentially unallocated or NULL array.\n");
      // Attempt to free what might exist, but be cautious
      if (array) {
          if (array[0]) { // Check if row pointers were allocated
             // We don't have a direct pointer to the contiguous data block
             // saved separately in this allocation scheme. The allocation relies
             // on array[0][0] pointing to it. If array[0] was freed first,
             // accessing array[0][0] is unsafe.
             // The allocation scheme where the contiguous data block is not
             // stored separately makes safe deallocation tricky if freeing
             // happens out of order or if parts are NULL.
 
             // A SAFER ALLOCATION/DEALLOCATION would store the data pointer separately.
             // Given the current allocation scheme, the order MUST be:
             // 1. Free the data block (pointed to by array[0][0])
             // 2. Free the row pointer block (pointed to by array[0])
             // 3. Free the layer pointer block (pointed to by array)
 
             // Let's assume the allocation was successful and pointers are valid.
             // Get pointer to the contiguous data block *before* freeing row pointers
             Cmpnts *dataContiguous = array[0][0];
             ierr = PetscFree(dataContiguous); CHKERRQ(ierr); // Free data block
 
             // Now free the row pointers block
             ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers
 
          }
          // Finally, free the array of layer pointers
          ierr = PetscFree(array); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0); // Return gracefully if input was NULL initially
  }
 
 
  // --- Standard Deallocation (assuming valid allocation) ---
 
  /* 1. Free the contiguous block of Cmpnts structures.
     The starting address was stored in array[0][0] by Allocate3DArrayVector. */
  ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
 
  /* 2. Free the contiguous block of row pointers.
     The starting address was stored in array[0]. */
  ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
 
  /* 3. Free the layer pointer array.
     The starting address is 'array' itself. */
  ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
 
  PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "GetOwnedCellRange"
/**
 * @brief Internal helper implementation: `GetOwnedCellRange()`.
 * @details Local to this translation unit.
 */
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes,
                                 PetscInt dim,
                                 PetscInt *xs_cell_global_out,
                                 PetscInt *xm_cell_local_out)
{
    PetscErrorCode ierr = 0; // Standard PETSc error code, not explicitly set here but good practice.
    PetscInt xs_node_global_rank;   // Global index of the first node owned by this rank in the specified dimension.
    PetscInt num_nodes_owned_rank;  // Number of nodes owned by this rank in this dimension (local count, excluding ghosts).
    PetscInt GlobalNodesInDim_from_info; // Total number of DA points in this dimension, from DMDALocalInfo.
 
    PetscFunctionBeginUser;
 
    // --- 1. Input Validation ---
    if (!info_nodes || !xs_cell_global_out || !xm_cell_local_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer passed to GetOwnedCellRange.");
    }
 
    // --- 2. Extract Node Ownership and Global Dimension Information from DMDALocalInfo ---
    if (dim == 0) { // I-direction
        xs_node_global_rank      = info_nodes->xs;
        num_nodes_owned_rank     = info_nodes->xm;
        GlobalNodesInDim_from_info = info_nodes->mx;
    } else if (dim == 1) { // J-direction
        xs_node_global_rank      = info_nodes->ys;
        num_nodes_owned_rank     = info_nodes->ym;
        GlobalNodesInDim_from_info = info_nodes->my;
    } else if (dim == 2) { // K-direction
        xs_node_global_rank      = info_nodes->zs;
        num_nodes_owned_rank     = info_nodes->zm;
        GlobalNodesInDim_from_info = info_nodes->mz;
    } else {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d in GetOwnedCellRange. Must be 0, 1, or 2.", dim);
    }
 
    // --- 3. Correct for User-Defined Ghost Node ---
    // Per the function's contract (@warning), the DA size includes an extra, non-physical
    // node. We subtract 1 to get the true number of physical nodes for cell calculations.
    const PetscInt physical_nodes_in_dim = GlobalNodesInDim_from_info - 1;
 
    // --- 4. Handle Edge Cases for Physical Domain Size ---
    // If the physical domain has 0 or 1 node, no cells can be formed.
    if (physical_nodes_in_dim <= 1) {
        *xs_cell_global_out = xs_node_global_rank; // Still report the rank's starting node
        *xm_cell_local_out  = 0;                    // But 0 cells
        PetscFunctionReturn(0);
    }
 
    // --- 5. Determine Cell Ownership Based on Corrected Node Ownership ---
    // The first cell this rank *could* define has its origin at the first node this rank owns.
    *xs_cell_global_out = xs_node_global_rank;
 
    // If the rank owns no nodes in this dimension, it can't form any cell origins.
    if (num_nodes_owned_rank == 0) {
        *xm_cell_local_out = 0;
    } else {
        // --- BUG FIX APPLIED HERE ---
        // The previous logic incorrectly assumed a cell's end node (N_{k+1}) must be on the
        // same rank as its origin node (N_k). The correct logic is to find the intersection
        // between the nodes this rank owns and the nodes that are valid origins globally.
 
        // The first node owned by the rank is its first potential origin.
        PetscInt first_owned_origin = xs_node_global_rank;
 
        // The absolute last node owned by this rank. Any node up to and including this one
        // is a potential cell origin from this rank's perspective.
        PetscInt last_node_owned_by_rank = xs_node_global_rank + num_nodes_owned_rank - 1;
 
        // The absolute last node in the entire PHYSICAL domain that can serve as a cell origin.
        // If there are `N` physical nodes (0 to N-1), this index is `N-2`.
        PetscInt last_possible_origin_global_idx = physical_nodes_in_dim - 2;
 
        // The actual last origin this rank can provide is the *minimum* of what it owns
        // and what is globally possible. This correctly handles both ranks in the middle of
        // the domain and the very last rank.
        PetscInt actual_last_origin_this_rank_can_form = PetscMin(last_node_owned_by_rank, last_possible_origin_global_idx);
 
        // If the first potential origin this rank owns is already beyond the actual last
        // origin it can form, then this rank forms no valid cell origins. This happens if
        // the rank only owns the very last physical node.
        if (first_owned_origin > actual_last_origin_this_rank_can_form) {
            *xm_cell_local_out = 0;
        } else {
            // The number of cells is the count of valid origins this rank owns.
            // (Count = Last Index - First Index + 1)
            *xm_cell_local_out = actual_last_origin_this_rank_can_form - first_owned_origin + 1;
        }
    }
 
    PetscFunctionReturn(ierr);
}
 
#undef __FUNCT__
#define __FUNCT__ "ComputeAndStoreNeighborRanks"
/**
 * @brief Internal helper implementation: `ComputeAndStoreNeighborRanks()`.
 * @details Local to this translation unit.
 */
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
{
    PetscErrorCode    ierr;
    PetscMPIInt       rank;
    PetscMPIInt       size; // MPI communicator size
    const PetscMPIInt *neighbor_ranks_ptr; // Pointer to raw neighbor data from PETSc
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); // Get MPI size for validation
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Computing DMDA neighbor ranks.\n", rank);
 
    if (!user || !user->da) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx or user->da is NULL in ComputeAndStoreNeighborRanks.");
    }
 
    // Get the neighbor information from the DMDA
    // neighbor_ranks_ptr will point to an internal PETSc array of 27 ranks.
    ierr = DMDAGetNeighbors(user->da, &neighbor_ranks_ptr); CHKERRQ(ierr);
 
    // Log the raw values from DMDAGetNeighbors for boundary-relevant directions for debugging
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "[Rank %d]Raw DMDAGetNeighbors: xm_raw=%d, xp_raw=%d, ym_raw=%d, yp_raw=%d, zm_raw=%d, zp_raw=%d. MPI_PROC_NULL is %d.\n",
                   rank,
                   neighbor_ranks_ptr[12], neighbor_ranks_ptr[14],
                   neighbor_ranks_ptr[10], neighbor_ranks_ptr[16],
                   neighbor_ranks_ptr[4],  neighbor_ranks_ptr[22],
                   (int)MPI_PROC_NULL);
 
    // PETSc standard indices for 3D face neighbors from the 27-point stencil:
    // Index = k_offset*9 + j_offset*3 + i_offset (where offsets -1,0,1 map to 0,1,2)
    // Center: (i_off=1, j_off=1, k_off=1) => 1*9 + 1*3 + 1 = 13
    // X-min:  (i_off=0, j_off=1, k_off=1) => 1*9 + 1*3 + 0 = 12
    // X-plus: (i_off=2, j_off=1, k_off=1) => 1*9 + 1*3 + 2 = 14
    // Y-min:  (i_off=1, j_off=0, k_off=1) => 1*9 + 0*3 + 1 = 10
    // Y-plus: (i_off=1, j_off=2, k_off=1) => 1*9 + 2*3 + 1 = 16
    // Z-min:  (i_off=1, j_off=1, k_off=0) => 0*9 + 1*3 + 1 = 4
    // Z-plus: (i_off=1, j_off=1, k_off=2) => 2*9 + 1*3 + 1 = 22
 
    if (neighbor_ranks_ptr[13] != rank) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Rank %d: DMDAGetNeighbors center index (13) is %d, expected current rank %d. Neighbor indexing might be non-standard or DMDA small.\n",
                  rank, neighbor_ranks_ptr[13], rank);
        // This warning is important. If the center isn't the current rank, the offsets are likely wrong.
        // However, PETSc should ensure this unless the DM is too small for a 3x3x3 stencil.
    }
 
    // Assign and sanitize each neighbor rank
    PetscMPIInt temp_neighbor;
 
    temp_neighbor = neighbor_ranks_ptr[12]; // xm
    if (temp_neighbor < 0  || temp_neighbor >= size) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid xm neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_xm = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_xm = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[14]; // xp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid xp neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_xp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_xp = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[10]; // ym
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid ym neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_ym = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_ym = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[16]; // yp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        // The log for index 16 was "zm" in your output, should be yp
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid yp neighbor (raw index 16) %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_yp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_yp = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[4]; // zm
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid zm neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_zm = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_zm = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[22]; // zp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid zp neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
        user->neighbors.rank_zp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_zp = temp_neighbor;
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "[Rank %d] Stored user->neighbors: xm=%d, xp=%d, ym=%d, yp=%d, zm=%d, zp=%d\n", rank,
              user->neighbors.rank_xm, user->neighbors.rank_xp,
              user->neighbors.rank_ym, user->neighbors.rank_yp,
              user->neighbors.rank_zm, user->neighbors.rank_zp);
    PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); // Ensure logs are flushed
 
    // Note: neighbor_ranks_ptr memory is managed by PETSc, do not free it.
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetDMDAProcLayout"
/**
 * @brief Internal helper implementation: `SetDMDAProcLayout()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetDMDAProcLayout(DM dm, UserCtx *user)
{
    PetscErrorCode ierr;
    PetscMPIInt    size, rank;
    PetscInt       px = PETSC_DECIDE, py = PETSC_DECIDE, pz = PETSC_DECIDE;
    PetscBool      px_set = PETSC_FALSE, py_set = PETSC_FALSE, pz_set = PETSC_FALSE;
    SimCtx         *simCtx = user->simCtx;
 
    // Set no.of processors in direction 1
    if(simCtx->da_procs_x) {
      px_set = PETSC_TRUE;
      px = simCtx->da_procs_x;
    }
    // Set no.of processors in direction 2
    if(simCtx->da_procs_y) {
      py_set = PETSC_TRUE;
      py = simCtx->da_procs_y;
    }
    // Set no.of processors in direction 1
    if(simCtx->da_procs_z) {
      pz_set = PETSC_TRUE;
      pz = simCtx->da_procs_z;
    }
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size); CHKERRQ(ierr);
    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Configuring DMDA processor layout for %d total processes.\n", rank, size);
 
    // --- Validate User Input (Optional but Recommended) ---
    // Check if specified processor counts multiply to the total MPI size
    if (px_set && py_set && pz_set) {
        if (px * py * pz != size) {
             SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP,
                     "Specified processor layout %d x %d x %d = %d does not match MPI size %d",
                     px, py, pz, px * py * pz, size);
        }
         LOG_ALLOW(GLOBAL, LOG_INFO, "Using specified processor layout: %d x %d x %d\n", px, py, pz);
    } else if (px_set || py_set || pz_set) {
         // If only some are set, PETSC_DECIDE will be used for others
         LOG_ALLOW(GLOBAL, LOG_INFO, "Using partially specified processor layout: %d x %d x %d (PETSC_DECIDE for unspecified)\n", px, py, pz);
    } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Using fully automatic processor layout (PETSC_DECIDE x PETSC_DECIDE x PETSC_DECIDE)\n");
    }
     // Additional checks: Ensure px, py, pz are positive if set
     if ((px_set && px <= 0) || (py_set && py <= 0) || (pz_set && pz <= 0)) {
         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Specified processor counts must be positive.");
     }
 
 
    // --- Apply the layout to the DMDA ---
    ierr = DMDASetNumProcs(dm, px, py, pz); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Rank %d: DMDASetNumProcs called with px=%d, py=%d, pz=%d.\n", rank, px, py, pz);
 
    // --- Store the values in UserCtx (Optional) ---
    // Note: If PETSC_DECIDE was used, PETSc calculates the actual values during DMSetUp.
    // We store the *requested* values here. To get the *actual* values used,
    // you would need to call DMDAGetInfo after DMSetUp.
    /*
    if (user) {
        user->procs_x = px;
        user->procs_y = py;
        user->procs_z = pz;
    }
    */
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetupDomainRankInfo"
/**
 * @brief Implementation of \ref SetupDomainRankInfo().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see SetupDomainRankInfo()
 */
PetscErrorCode SetupDomainRankInfo(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscInt       nblk = simCtx->block_number;
    PetscInt       size = simCtx->size;
    BoundingBox    *final_bboxlist = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting full rank communication setup for %d block(s).\n", nblk);
 
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    // --- Step 1: Compute neighbor ranks (unchanged) ---
    for (int bi = 0; bi < nblk; bi++) {
        ierr = ComputeAndStoreNeighborRanks(&user_finest[bi]); CHKERRQ(ierr);
    }
    LOG_ALLOW(GLOBAL, LOG_INFO, "Neighbor ranks computed and stored for all blocks.\n");
 
    // --- Step 2: Allocate the final, unified list on ALL ranks ---
    // Every rank will build this list in parallel.
    ierr = PetscMalloc1(size * nblk, &final_bboxlist); CHKERRQ(ierr);
 
    // --- Step 3: Loop through each block, gather then broadcast its bbox list ---
    for (int bi = 0; bi < nblk; bi++) {
        // This is a temporary pointer for the current block's list.
        BoundingBox *block_bboxlist = NULL;
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Processing bounding boxes for block %d...\n", bi);
 
        // A) GATHER: On rank 0, block_bboxlist is allocated and filled. On others, it's NULL.
        ierr = GatherAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Gather complete for block %d.\n", bi);
 
        // B) BROADCAST: On non-root ranks, block_bboxlist is allocated. Then, the data
        //    from rank 0 is broadcast to all ranks. After this call, ALL ranks have
        //    an identical, complete copy of the bounding boxes for the current block.
        ierr = BroadcastAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Broadcast complete for block %d.\n", bi);
 
        // C) ASSEMBLE: Every rank now copies the data for this block into the
        //    correct segment of its final, unified list.
        for (int r = 0; r < size; r++) {
            // The layout is [r0b0, r1b0, ..., r(size-1)b0, r0b1, r1b1, ...]
            final_bboxlist[bi * size + r] = block_bboxlist[r];
        }
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Assembly into final list complete for block %d.\n", bi);
 
        // D) CLEANUP: Free the temporary list for this block on ALL ranks before the next iteration.
        //    Your helper functions use malloc, so we must use free.
        free(block_bboxlist);
    }
 
    // --- Step 4: Assign the final pointer and run the last setup step ---
    simCtx->bboxlist = final_bboxlist;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Final unified bboxlist created on all ranks and stored in SimCtx.\n");
 
    ierr = SetupDomainCellDecompositionMap(&user_finest[0]); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Domain Cell Composition set and broadcasted.\n");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "SetupDomainRankInfo: Completed successfully.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "Contra2Cart"
/**
 * @brief Internal helper implementation: `Contra2Cart()`.
 * @details Local to this translation unit.
 */
PetscErrorCode Contra2Cart(UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Cmpnts       ***lcsi_arr, ***leta_arr, ***lzet_arr; // Local metric arrays
    Cmpnts       ***lucont_arr;                       // Local contravariant velocity array
    Cmpnts       ***gucat_arr;                        // Global Cartesian velocity array
    PetscReal    ***lnvert_arr;                       // Local Nvert array
    PetscReal    ***laj_arr;                          // Local Jacobian Determinant inverse array
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting Contravariant-to-Cartesian velocity transformation.\n");
 
    // --- 1. Get DMDA Info and Check for Valid Inputs ---
    // All inputs (lUcont, lCsi, etc.) and outputs (Ucat) are on DMs from the UserCtx.
    // We get local info from fda, which governs the layout of most arrays here.
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    if (!user->lUcont || !user->lCsi || !user->lEta || !user->lZet || !user->lNvert || !user->Ucat) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Contra2Cart requires lUcont, lCsi/Eta/Zet, lNvert, and Ucat to be non-NULL.");
    }
 
 
    // --- 2. Get Read-Only Array Access to Local Input Vectors (with ghosts) ---
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi,   &lcsi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta,   &leta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet,   &lzet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lNvert, &lnvert_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj, &laj_arr); CHKERRQ(ierr);
 
    // --- 3. Get Write-Only Array Access to the Global Output Vector ---
    // We compute for local owned cells and write into the global vector.
    // PETSc handles mapping the global indices to the correct local memory locations.
    ierr = DMDAVecGetArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
 
 
    // --- 4. Define Loop Bounds for INTERIOR Cells ---
    // We use adjusted bounds to avoid calculating Ucat on the physical domain boundaries,
    // as these are typically set explicitly by boundary condition functions.
    // The stencils use indices like i-1, j-1, k-1, so we must start loops at least at index 1.
    PetscInt i_start = (info.xs == 0) ? info.xs + 1 : info.xs;
    PetscInt i_end   = (info.xs + info.xm == info.mx) ? info.xs + info.xm - 1 : info.xs + info.xm;
 
    PetscInt j_start = (info.ys == 0) ? info.ys + 1 : info.ys;
    PetscInt j_end   = (info.ys + info.ym == info.my) ? info.ys + info.ym - 1 : info.ys + info.ym;
 
    PetscInt k_start = (info.zs == 0) ? info.zs + 1 : info.zs;
    PetscInt k_end   = (info.zs + info.zm == info.mz) ? info.zs + info.zm - 1 : info.zs + info.zm;
 
    // --- 5. Main Computation Loop ---
    // Loops over the GLOBAL indices of interior cells owned by this rank.
    for (PetscInt k_cell = k_start; k_cell < k_end; ++k_cell) {
        for (PetscInt j_cell = j_start; j_cell < j_end; ++j_cell) {
            for (PetscInt i_cell = i_start; i_cell < i_end; ++i_cell) {
 
                // Check if the cell is a fluid cell (not solid/blanked)
          //    if (lnvert_arr[k_cell][j_cell][i_cell] > 0.1) continue; // Skip solid/blanked cells
 
                // Transformation matrix [mat] is the metric tensor at the cell center,
                // estimated by averaging metrics from adjacent faces.
                PetscReal mat[3][3];
 
        //  PetscReal aj_center = laj_arr[k_cell+1][j_cell+1][i_cell+1];
        
                mat[0][0] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].x + lcsi_arr[k_cell][j_cell][i_cell].x); //* aj_center;
                mat[0][1] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].y + lcsi_arr[k_cell][j_cell][i_cell].y); //* aj_center;
                mat[0][2] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].z + lcsi_arr[k_cell][j_cell][i_cell].z); //* aj_center;
 
                mat[1][0] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].x + leta_arr[k_cell][j_cell][i_cell].x); //* aj_center;
                mat[1][1] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].y + leta_arr[k_cell][j_cell][i_cell].y); //* aj_center;
                mat[1][2] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].z + leta_arr[k_cell][j_cell][i_cell].z); //* aj_center;
 
                mat[2][0] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].x + lzet_arr[k_cell][j_cell][i_cell].x); //* aj_center;
                mat[2][1] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].y + lzet_arr[k_cell][j_cell][i_cell].y); //* aj_center;
                mat[2][2] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].z + lzet_arr[k_cell][j_cell][i_cell].z); //* aj_center;
    
                // Contravariant velocity vector `q` at the cell center,
                // estimated by averaging face-based contravariant velocities.
                PetscReal q[3];
                q[0] = 0.5 * (lucont_arr[k_cell][j_cell][i_cell-1].x + lucont_arr[k_cell][j_cell][i_cell].x); // U¹ at cell center
                q[1] = 0.5 * (lucont_arr[k_cell][j_cell-1][i_cell].y + lucont_arr[k_cell][j_cell][i_cell].y); // U² at cell center
                q[2] = 0.5 * (lucont_arr[k_cell-1][j_cell][i_cell].z + lucont_arr[k_cell][j_cell][i_cell].z); // U³ at cell center
 
                // Solve the 3x3 system `mat * ucat = q` using Cramer's rule.
                PetscReal det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
                                mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
                                mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
 
          if (PetscAbsReal(det) < 1.0e-18) {
              SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Transformation matrix determinant is near zero at cell (%d,%d,%d) \n", i_cell, j_cell, k_cell);
                        }
 
                PetscReal det_inv = 1.0 / det;
 
                PetscReal det0 = q[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
                                 q[1] * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) +
                                 q[2] * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
 
                PetscReal det1 = -q[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
                                  q[1] * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) -
                                  q[2] * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
 
                PetscReal det2 = q[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) -
                                 q[1] * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) +
                                 q[2] * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
 
                // Store computed Cartesian velocity in the GLOBAL Ucat array at the
                // array index corresponding to the cell's origin node.
                gucat_arr[k_cell][j_cell][i_cell].x = det0 * det_inv;
                gucat_arr[k_cell][j_cell][i_cell].y = det1 * det_inv;
                gucat_arr[k_cell][j_cell][i_cell].z = det2 * det_inv;
            }
        }
    }
 
    // --- 6. Restore Array Access ---
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi,   &lcsi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta,   &leta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet,   &lzet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lNvert, &lnvert_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj, &laj_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Completed Contravariant-to-Cartesian velocity transformation. \n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "Cart2Contra"
/**
 * @brief Convert a spatially varying Cartesian velocity field to contravariant fluxes.
 */
PetscErrorCode Cart2Contra(UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo info;
    const Cmpnts ***ucat_arr, ***csi_arr, ***eta_arr, ***zet_arr;
    Cmpnts ***ucont_arr;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
 
    const PetscInt i_start = PetscMax(info.xs, 1);
    const PetscInt j_start = PetscMax(info.ys, 1);
    const PetscInt k_start = PetscMax(info.zs, 1);
    const PetscInt i_end = PetscMin(info.xs + info.xm, info.mx - 1);
    const PetscInt j_end = PetscMin(info.ys + info.ym, info.my - 1);
    const PetscInt k_end = PetscMin(info.zs + info.zm, info.mz - 1);
 
    for (PetscInt k = k_start; k < k_end; k++) {
        for (PetscInt j = j_start; j < j_end; j++) {
            for (PetscInt i = i_start; i < i_end; i++) {
                const Cmpnts u_xi = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j][i + 1].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j][i + 1].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j][i + 1].z)
                };
                const Cmpnts u_eta = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j + 1][i].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j + 1][i].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j + 1][i].z)
                };
                const Cmpnts u_zeta = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k + 1][j][i].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k + 1][j][i].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k + 1][j][i].z)
                };
                ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u_xi.x + csi_arr[k][j][i].y * u_xi.y + csi_arr[k][j][i].z * u_xi.z;
                ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u_eta.x + eta_arr[k][j][i].y * u_eta.y + eta_arr[k][j][i].z * u_eta.z;
                ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u_zeta.x + zet_arr[k][j][i].y * u_zeta.y + zet_arr[k][j][i].z * u_zeta.z;
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "UniformCart2Contra"
/*
 * Converts a uniform Cartesian velocity (u,v,w) to contravariant fluxes in Ucont.
 * Computes the dot product of the physical velocity with each face-area vector:
 *   U^xi  = csi · (u,v,w),  U^eta = eta · (u,v,w),  U^zeta = zet · (u,v,w).
 * Writes to all owned nodes (xs..xe, ys..ye, zs..ze); boundary ghosts are
 * overwritten later by ApplyBoundaryConditions.
 * The public header owns the rendered API contract.
 */
PetscErrorCode UniformCart2Contra(UserCtx *user, PetscReal u, PetscReal v, PetscReal w)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    DMDALocalInfo    info;
    Cmpnts         ***ucont_arr;
    const Cmpnts   ***csi_arr, ***eta_arr, ***zet_arr;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr);        CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr);       CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr);       CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr);       CHKERRQ(ierr);
 
    const PetscInt xs = info.xs, xe = info.xs + info.xm;
    const PetscInt ys = info.ys, ye = info.ys + info.ym;
    const PetscInt zs = info.zs, ze = info.zs + info.zm;
 
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u + csi_arr[k][j][i].y * v + csi_arr[k][j][i].z * w;
                ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u + eta_arr[k][j][i].y * v + eta_arr[k][j][i].z * w;
                ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u + zet_arr[k][j][i].y * v + zet_arr[k][j][i].z * w;
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr);    CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Cart2Contra: set Ucont from Cartesian (%.3f, %.3f, %.3f).\n",
              (double)u, (double)v, (double)w);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetupDomainCellDecompositionMap"
/**
 * @brief Internal helper implementation: `SetupDomainCellDecompositionMap()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  local_node_info;
    RankCellInfo   my_cell_info;
    PetscMPIInt    rank, size;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation and MPI Info ---
    if (!user) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL in SetupDomainCellDecompositionMap.");
    }
    if (!user->da) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->da is not initialized in SetupDomainCellDecompositionMap.");
    }
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Setting up domain cell decomposition map for %d ranks.\n", size);
 
    // --- 2. Determine Local Cell Ownership ---
    // Get the local node ownership information from the primary DMDA.
    ierr = DMDAGetLocalInfo(user->da, &local_node_info); CHKERRQ(ierr);
 
    // Use the robust helper function to convert node ownership to cell ownership.
    // A cell's index is defined by its origin node.
     
    ierr = GetOwnedCellRange(&local_node_info, 0, &my_cell_info.xs_cell, &my_cell_info.xm_cell); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&local_node_info, 1, &my_cell_info.ys_cell, &my_cell_info.ym_cell); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&local_node_info, 2, &my_cell_info.zs_cell, &my_cell_info.zm_cell); CHKERRQ(ierr);
    
    // Log the calculated local ownership for debugging purposes.
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Owns cells: i[%d, %d), j[%d, %d), k[%d, %d)\n",
              rank, my_cell_info.xs_cell, my_cell_info.xs_cell + my_cell_info.xm_cell,
              my_cell_info.ys_cell, my_cell_info.ys_cell + my_cell_info.ym_cell,
              my_cell_info.zs_cell, my_cell_info.zs_cell + my_cell_info.zm_cell);
 
    // --- 3. Allocate and Distribute the Global Map ---
    // Allocate memory for the global map that will hold information from all ranks.
    ierr = PetscMalloc1(size, &user->RankCellInfoMap); CHKERRQ(ierr);
 
    // Perform the collective communication to gather the `RankCellInfo` struct from every rank.
    // Each rank sends its `my_cell_info` and receives the complete array in `user->RankCellInfoMap`.
    // We use MPI_BYTE to ensure portability across different systems and struct padding.
    ierr = MPI_Allgather(&my_cell_info, sizeof(RankCellInfo), MPI_BYTE,
                         user->RankCellInfoMap, sizeof(RankCellInfo), MPI_BYTE,
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Domain cell decomposition map created and distributed successfully.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "BinarySearchInt64"
/**
 * @brief Implementation of \ref BinarySearchInt64().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see BinarySearchInt64()
 */
PetscErrorCode BinarySearchInt64(PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
{
    PetscInt low = 0, high = n - 1;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation ---
    if (!found) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer 'found' is NULL in PetscBinarySearchInt64.");
    }
    if (n > 0 && !arr) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input array 'arr' is NULL for n > 0.");
    }
    
    // Initialize output
    *found = PETSC_FALSE;
 
    // --- 2. Binary Search Algorithm ---
    while (low <= high) {
        // Use this form to prevent potential integer overflow on very large arrays
        PetscInt mid = low + (high - low) / 2;
 
        if (arr[mid] == key) {
            *found = PETSC_TRUE; // Key found!
            break;               // Exit the loop
        }
        
        if (arr[mid] < key) {
            low = mid + 1; // Search in the right half
        } else {
            high = mid - 1; // Search in the left half
        }
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
 
/**
 * @brief Convert logical indices into the flattened global index used by setup helpers.
 */
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
{
  PetscInt nidx;
  DMDALocalInfo info = user->info;
 
  PetscInt  mx = info.mx, my = info.my;
  
  AO ao;
  DMDAGetAO(user->da, &ao);
  nidx=i+j*mx+k*mx*my;
  
  AOApplicationToPetsc(ao,1,&nidx);
  
  return (nidx);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "ComputeDivergence"
/**
 * @brief Implementation of \ref ComputeDivergence().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see ComputeDivergence()
 */
 
PetscErrorCode ComputeDivergence(UserCtx *user )
{
  DM        da = user->da, fda = user->fda;
  DMDALocalInfo info = user->info;
 
  PetscInt ti = user->simCtx->step;
  
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt  lxs, lys, lzs, lxe, lye, lze;
  PetscInt  i, j, k;
 
  Vec       Div;
  PetscReal ***div, ***aj, ***nvert,***p;
  Cmpnts    ***ucont;
  PetscReal maxdiv;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
 
  DMDAVecGetArray(fda,user->lUcont, &ucont);
  DMDAVecGetArray(da, user->lAj, &aj);
  VecDuplicate(user->P, &Div);
  DMDAVecGetArray(da, Div, &div);
  DMDAVecGetArray(da, user->lNvert, &nvert);
   DMDAVecGetArray(da, user->P, &p);
  for (k=lzs; k<lze; k++) {
         for (j=lys; j<lye; j++){
      for (i=lxs; i<lxe; i++) {
        if (k==10 && j==10 && i==1){
          LOG_ALLOW(LOCAL,LOG_INFO,"Pressure[10][10][1] =  %f | Pressure[10][10][0] =  %f  \n ",p[k][j][i],p[k][j][i-1]);
              }
 
        if (k==10 && j==10 && i==mx-3)
              LOG_ALLOW(LOCAL,LOG_INFO,"Pressure[10][10][%d] =  %f | Pressure[10][10][%d] =  %f  \n ",mx-2,p[k][j][mx-2],mx-1,p[k][j][mx-1]);
        }
    }
    }
   DMDAVecRestoreArray(da, user->P, &p);
 
 
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    maxdiv = fabs((ucont[k][j][i].x - ucont[k][j][i-1].x +
               ucont[k][j][i].y - ucont[k][j-1][i].y +
               ucont[k][j][i].z - ucont[k-1][j][i].z)*aj[k][j][i]);
    if (nvert[k][j][i] + nvert[k+1][j][i] + nvert[k-1][j][i] +
        nvert[k][j+1][i] + nvert[k][j-1][i] +
        nvert[k][j][i+1] + nvert[k][j][i-1] > 0.1) maxdiv = 0.;
    div[k][j][i] = maxdiv;
 
      }
    }
  }
 
  if (zs==0) {
    k=0;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (ze == mz) {
    k=mz-1;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (xs==0) {
    i=0;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (xe==mx) {
    i=mx-1;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    div[k][j][i] = 0;
      }
    }
  }
 
  if (ys==0) {
    j=0;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (ye==my) {
    j=my-1;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
  DMDAVecRestoreArray(da, Div, &div);
  PetscInt MaxFlatIndex;
  
  VecMax(Div, &MaxFlatIndex, &maxdiv);
 
  LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d]] The Maximum Divergence is %e at flat index %d.\n",ti,maxdiv,MaxFlatIndex); 
 
  user->simCtx->MaxDivFlatArg = MaxFlatIndex;
  user->simCtx->MaxDiv = maxdiv;
  
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (Gidx(i,j,k,user) == MaxFlatIndex) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d] The Maximum Divergence(%e) is at location [%d][%d][%d]. \n", ti, maxdiv,k,j,i);
      user->simCtx->MaxDivz = k;
      user->simCtx->MaxDivy = j;
      user->simCtx->MaxDivx = i;
    }
      }
    }
  }
 
  
 DMDAVecRestoreArray(da, user->lNvert, &nvert);
 DMDAVecRestoreArray(fda, user->lUcont, &ucont);
 DMDAVecRestoreArray(da, user->lAj, &aj);
 VecDestroy(&Div);
 
 PROFILE_FUNCTION_END;
 PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InitializeRandomGenerators"
 
/**
 * @brief Implementation of \ref InitializeRandomGenerators().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see InitializeRandomGenerators()
 */
PetscErrorCode InitializeRandomGenerators(UserCtx* user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz) {
    PetscErrorCode ierr;  // Error code for PETSc functions
    PetscMPIInt rank;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
    // Initialize RNG for x-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randx); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randx), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randx, user->bbox.min_coords.x, user->bbox.max_coords.x); CHKERRQ(ierr);
    /* Offsets reproduce the historical fixed seeds (12345, 67890, 54321 plus rank) at the
       default base seed, so existing placements are unchanged. */
    const unsigned long base_seed = (unsigned long)user->simCtx->particleRandomSeed;
    ierr = PetscRandomSetSeed(*randx, base_seed + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randx); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for X-axis.\n",rank);
 
    // Initialize RNG for y-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randy); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randy), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randy, user->bbox.min_coords.y, user->bbox.max_coords.y); CHKERRQ(ierr);
    ierr = PetscRandomSetSeed(*randy, base_seed + 55545UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randy); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Y-axis.\n",rank);
 
    // Initialize RNG for z-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randz); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randz), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randz, user->bbox.min_coords.z, user->bbox.max_coords.z); CHKERRQ(ierr);
    ierr = PetscRandomSetSeed(*randz, base_seed + 41976UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randz); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Z-axis.\n",rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InitializeLogicalSpaceRNGs"
/**
 * @brief Internal helper implementation: `InitializeLogicalSpaceRNGs()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InitializeLogicalSpaceRNGs(PetscInt base_seed, PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k) {
    PetscErrorCode ierr;
    PetscMPIInt rank;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- RNG for i-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_i); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_i), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_i, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    /* Offsets reproduce the historical seeds (202401..202403 plus rank) at the default base. */
    ierr = PetscRandomSetSeed(*rand_logic_i, (unsigned long)base_seed + 190056UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_i); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for i-logical dimension [0,1).\n",rank);
 
    // --- RNG for j-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_j); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_j), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_j, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    ierr = PetscRandomSetSeed(*rand_logic_j, (unsigned long)base_seed + 190057UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_j); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for j-logical dimension [0,1).\n",rank);
 
    // --- RNG for k-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_k); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_k), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_k, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    ierr = PetscRandomSetSeed(*rand_logic_k, (unsigned long)base_seed + 190058UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_k); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for k-logical dimension [0,1).\n",rank);
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InitializeBrownianRNG"
/**
 * @brief Internal helper implementation: `InitializeBrownianRNG()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InitializeBrownianRNG(SimCtx *simCtx) {
    PetscErrorCode ierr;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // 1. Create the generator (stored in SimCtx, not UserCtx, as it is global physics)
    ierr = PetscRandomCreate(PETSC_COMM_WORLD, &simCtx->BrownianMotionRNG); CHKERRQ(ierr);
    ierr = PetscRandomSetType(simCtx->BrownianMotionRNG, PETSCRAND48); CHKERRQ(ierr);
 
    // 2. CRITICAL: Set interval to [0, 1). 
    // This is required for the Gaussian math to work.
    ierr = PetscRandomSetInterval(simCtx->BrownianMotionRNG, 0.0, 1.0); CHKERRQ(ierr);
 
    // 3. Seed from the configured base, the rank, and the start step. It was seeded from
    //    the wall clock, so no two runs of the same inputs agreed. The start step keeps a
    //    restart from replaying the increments the first segment already drew, which would
    //    correlate them; the same restart is still reproducible.
    unsigned long seed = (unsigned long)simCtx->particleRandomSeed
                       + (unsigned long)rank * 987654321UL
                       + (unsigned long)simCtx->StartStep * 1000003UL;
    ierr = PetscRandomSetSeed(simCtx->BrownianMotionRNG, seed); CHKERRQ(ierr);
    ierr = PetscRandomSeed(simCtx->BrownianMotionRNG); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE, "[Rank %d] Initialized Brownian Physics RNG.\n", rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
/////////////// DERIVATIVE CALCULATION HELPERS ///////////////
 
#undef __FUNCT__
#define __FUNCT__ "TransformScalarDerivativesToPhysical"
/**
 * @brief Implementation of \ref TransformScalarDerivativesToPhysical().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see TransformScalarDerivativesToPhysical()
 */
 void TransformScalarDerivativesToPhysical(PetscReal jacobian, 
                                                 Cmpnts csi_metrics, 
                                                 Cmpnts eta_metrics, 
                                                 Cmpnts zet_metrics,
                                                 PetscReal dPhi_dcsi, 
                                                 PetscReal dPhi_deta, 
                                                 PetscReal dPhi_dzet,
                                                 Cmpnts *gradPhi)
{
    // Gradient X component
    gradPhi->x = jacobian * (dPhi_dcsi * csi_metrics.x + dPhi_deta * eta_metrics.x + dPhi_dzet * zet_metrics.x);
    
    // Gradient Y component
    gradPhi->y = jacobian * (dPhi_dcsi * csi_metrics.y + dPhi_deta * eta_metrics.y + dPhi_dzet * zet_metrics.y);
    
    // Gradient Z component
    gradPhi->z = jacobian * (dPhi_dcsi * csi_metrics.z + dPhi_deta * eta_metrics.z + dPhi_dzet * zet_metrics.z);
}
 
#undef __FUNCT__
#define __FUNCT__ "TransformDerivativesToPhysical"
/**
 * @brief Transform contravariant vector derivatives into physical Cartesian derivatives.
 */
static void TransformDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics,
                                           Cmpnts deriv_csi, Cmpnts deriv_eta, Cmpnts deriv_zet,
                                           Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
{
    // Derivatives of the first component (u)
    dudx->x = jacobian * (deriv_csi.x * csi_metrics.x + deriv_eta.x * eta_metrics.x + deriv_zet.x * zet_metrics.x);
    dudx->y = jacobian * (deriv_csi.x * csi_metrics.y + deriv_eta.x * eta_metrics.y + deriv_zet.x * zet_metrics.y);
    dudx->z = jacobian * (deriv_csi.x * csi_metrics.z + deriv_eta.x * eta_metrics.z + deriv_zet.x * zet_metrics.z);
    // Derivatives of the second component (v)
    dvdx->x = jacobian * (deriv_csi.y * csi_metrics.x + deriv_eta.y * eta_metrics.x + deriv_zet.y * zet_metrics.x);
    dvdx->y = jacobian * (deriv_csi.y * csi_metrics.y + deriv_eta.y * eta_metrics.y + deriv_zet.y * zet_metrics.y);
    dvdx->z = jacobian * (deriv_csi.y * csi_metrics.z + deriv_eta.y * eta_metrics.z + deriv_zet.y * zet_metrics.z);
    // Derivatives of the third component (w)
    dwdx->x = jacobian * (deriv_csi.z * csi_metrics.x + deriv_eta.z * eta_metrics.x + deriv_zet.z * zet_metrics.x);
    dwdx->y = jacobian * (deriv_csi.z * csi_metrics.y + deriv_eta.z * eta_metrics.y + deriv_zet.z * zet_metrics.y);
    dwdx->z = jacobian * (deriv_csi.z * csi_metrics.z + deriv_eta.z * eta_metrics.z + deriv_zet.z * zet_metrics.z);
}
 
#undef __FUNCT__
#define __FUNCT__ "ComputeScalarFieldDerivatives"
/**
 * @brief Internal helper implementation: `ComputeScalarFieldDerivatives()`.
 * @details Local to this translation unit.
 */
PetscErrorCode ComputeScalarFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, 
                                             PetscReal ***field_data, Cmpnts *grad)
{
    PetscErrorCode ierr;
    Cmpnts    ***csi, ***eta, ***zet;
    PetscReal ***jac;
    PetscReal d_csi, d_eta, d_zet;
 
    PetscFunctionBeginUser;
 
    // 1. Get read-only access to metrics
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    // 2. Compute derivatives in computational space (Central Difference)
    //    Assumes ghosts are available at i+/-1
    d_csi = 0.5 * (field_data[k][j][i+1] - field_data[k][j][i-1]);
    d_eta = 0.5 * (field_data[k][j+1][i] - field_data[k][j-1][i]);
    d_zet = 0.5 * (field_data[k+1][j][i] - field_data[k-1][j][i]);
 
    // 3. Transform to physical space
    TransformScalarDerivativesToPhysical(jac[k][j][i], 
                                         csi[k][j][i], eta[k][j][i], zet[k][j][i],
                                         d_csi, d_eta, d_zet,
                                         grad);
 
    // 4. Restore arrays
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ComputeVectorFieldDerivatives"
/**
 * @brief Internal helper implementation: `ComputeVectorFieldDerivatives()`.
 * @details Local to this translation unit.
 */
PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data,
                                             Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
{
    PetscErrorCode ierr;
    Cmpnts ***csi, ***eta, ***zet;
    PetscReal ***jac;
    PetscFunctionBeginUser;
 
    // 1. Get read-only access to the necessary metric data arrays
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    // 2. Calculate derivatives in computational space using central differencing
    Cmpnts deriv_csi, deriv_eta, deriv_zet;
    deriv_csi.x = (field_data[k][j][i+1].x - field_data[k][j][i-1].x) * 0.5;
    deriv_csi.y = (field_data[k][j][i+1].y - field_data[k][j][i-1].y) * 0.5;
    deriv_csi.z = (field_data[k][j][i+1].z - field_data[k][j][i-1].z) * 0.5;
 
    deriv_eta.x = (field_data[k][j+1][i].x - field_data[k][j-1][i].x) * 0.5;
    deriv_eta.y = (field_data[k][j+1][i].y - field_data[k][j-1][i].y) * 0.5;
    deriv_eta.z = (field_data[k][j+1][i].z - field_data[k][j-1][i].z) * 0.5;
 
    deriv_zet.x = (field_data[k+1][j][i].x - field_data[k-1][j][i].x) * 0.5;
    deriv_zet.y = (field_data[k+1][j][i].y - field_data[k-1][j][i].y) * 0.5;
    deriv_zet.z = (field_data[k+1][j][i].z - field_data[k-1][j][i].z) * 0.5;
 
    // 3. Transform derivatives to physical space
    TransformDerivativesToPhysical(jac[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i],
                                   deriv_csi, deriv_eta, deriv_zet,
                                   dudx, dvdx, dwdx);
 
    // 4. Restore access to the PETSc data arrays
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}
 
//================================================================================
//
//                         MEMORY CLEANUP FUNCTIONS
//
//================================================================================
 
#undef __FUNCT__
#define __FUNCT__ "DestroyUserVectors"
/**
 * @brief Internal helper implementation: `DestroyUserVectors()`.
 * @details Local to this translation unit.
 */
PetscErrorCode DestroyUserVectors(UserCtx *user)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    // --- Group A: Primary Flow Fields (Always allocated at all levels) ---
    if (user->Ucont) { ierr = VecDestroy(&user->Ucont); CHKERRQ(ierr); }
    if (user->lUcont) { ierr = VecDestroy(&user->lUcont); CHKERRQ(ierr); }
    if (user->Ucat) { ierr = VecDestroy(&user->Ucat); CHKERRQ(ierr); }
    if (user->lUcat) { ierr = VecDestroy(&user->lUcat); CHKERRQ(ierr); }
    if (user->P) { ierr = VecDestroy(&user->P); CHKERRQ(ierr); }
    if (user->lP) { ierr = VecDestroy(&user->lP); CHKERRQ(ierr); }
    if (user->Nvert) { ierr = VecDestroy(&user->Nvert); CHKERRQ(ierr); }
    if (user->lNvert) { ierr = VecDestroy(&user->lNvert); CHKERRQ(ierr); }
 
    // --- Group A2: Derived Flow Fields (Conditional) ---
    if(user->Diffusivity) {ierr = VecDestroy(&user->Diffusivity); CHKERRQ(ierr);}
    if(user->lDiffusivity){ierr = VecDestroy(&user->lDiffusivity); CHKERRQ(ierr);}
    if(user->DiffusivityGradient){ierr = VecDestroy(&user->DiffusivityGradient); CHKERRQ(ierr);}
    if(user->lDiffusivityGradient){ierr = VecDestroy(&user->lDiffusivityGradient); CHKERRQ(ierr);}
    
    // --- Group B: Solver Work Vectors (All levels) ---
    if (user->Phi) { ierr = VecDestroy(&user->Phi); CHKERRQ(ierr); }
    if (user->lPhi) { ierr = VecDestroy(&user->lPhi); CHKERRQ(ierr); }
 
    // --- Group C: Time-Stepping Vectors (Finest level only) ---
    if (user->Ucont_o) { ierr = VecDestroy(&user->Ucont_o); CHKERRQ(ierr); }
    if (user->Ucont_rm1) { ierr = VecDestroy(&user->Ucont_rm1); CHKERRQ(ierr); }
    if (user->Ucat_o) { ierr = VecDestroy(&user->Ucat_o); CHKERRQ(ierr); }
    if (user->P_o) { ierr = VecDestroy(&user->P_o); CHKERRQ(ierr); }
    if (user->Nvert_o) { ierr = VecDestroy(&user->Nvert_o); CHKERRQ(ierr); }
    if (user->lUcont_o) { ierr = VecDestroy(&user->lUcont_o); CHKERRQ(ierr); }
    if (user->lUcont_rm1) { ierr = VecDestroy(&user->lUcont_rm1); CHKERRQ(ierr); }
    if (user->lNvert_o) { ierr = VecDestroy(&user->lNvert_o); CHKERRQ(ierr); }
 
    // --- Group D: Grid Metrics - Face Centered (All levels) ---
    if (user->Csi) { ierr = VecDestroy(&user->Csi); CHKERRQ(ierr); }
    if (user->Eta) { ierr = VecDestroy(&user->Eta); CHKERRQ(ierr); }
    if (user->Zet) { ierr = VecDestroy(&user->Zet); CHKERRQ(ierr); }
    if (user->Aj) { ierr = VecDestroy(&user->Aj); CHKERRQ(ierr); }
    if (user->lCsi) { ierr = VecDestroy(&user->lCsi); CHKERRQ(ierr); }
    if (user->lEta) { ierr = VecDestroy(&user->lEta); CHKERRQ(ierr); }
    if (user->lZet) { ierr = VecDestroy(&user->lZet); CHKERRQ(ierr); }
    if (user->lAj) { ierr = VecDestroy(&user->lAj); CHKERRQ(ierr); }
 
    // --- Group E: Grid Metrics - Face Centered (All levels) ---
    if (user->ICsi) { ierr = VecDestroy(&user->ICsi); CHKERRQ(ierr); }
    if (user->IEta) { ierr = VecDestroy(&user->IEta); CHKERRQ(ierr); }
    if (user->IZet) { ierr = VecDestroy(&user->IZet); CHKERRQ(ierr); }
    if (user->JCsi) { ierr = VecDestroy(&user->JCsi); CHKERRQ(ierr); }
    if (user->JEta) { ierr = VecDestroy(&user->JEta); CHKERRQ(ierr); }
    if (user->JZet) { ierr = VecDestroy(&user->JZet); CHKERRQ(ierr); }
    if (user->KCsi) { ierr = VecDestroy(&user->KCsi); CHKERRQ(ierr); }
    if (user->KEta) { ierr = VecDestroy(&user->KEta); CHKERRQ(ierr); }
    if (user->KZet) { ierr = VecDestroy(&user->KZet); CHKERRQ(ierr); }
    if (user->IAj) { ierr = VecDestroy(&user->IAj); CHKERRQ(ierr); }
    if (user->JAj) { ierr = VecDestroy(&user->JAj); CHKERRQ(ierr); }
    if (user->KAj) { ierr = VecDestroy(&user->KAj); CHKERRQ(ierr); }
    if (user->lICsi) { ierr = VecDestroy(&user->lICsi); CHKERRQ(ierr); }
    if (user->lIEta) { ierr = VecDestroy(&user->lIEta); CHKERRQ(ierr); }
    if (user->lIZet) { ierr = VecDestroy(&user->lIZet); CHKERRQ(ierr); }
    if (user->lJCsi) { ierr = VecDestroy(&user->lJCsi); CHKERRQ(ierr); }
    if (user->lJEta) { ierr = VecDestroy(&user->lJEta); CHKERRQ(ierr); }
    if (user->lJZet) { ierr = VecDestroy(&user->lJZet); CHKERRQ(ierr); }
    if (user->lKCsi) { ierr = VecDestroy(&user->lKCsi); CHKERRQ(ierr); }
    if (user->lKEta) { ierr = VecDestroy(&user->lKEta); CHKERRQ(ierr); }
    if (user->lKZet) { ierr = VecDestroy(&user->lKZet); CHKERRQ(ierr); }
    if (user->lIAj) { ierr = VecDestroy(&user->lIAj); CHKERRQ(ierr); }
    if (user->lJAj) { ierr = VecDestroy(&user->lJAj); CHKERRQ(ierr); }
    if (user->lKAj) { ierr = VecDestroy(&user->lKAj); CHKERRQ(ierr); }
 
    // --- Group F: Cell/Face Coordinates and Grid Spacing (All levels) ---
    if (user->Cent) { ierr = VecDestroy(&user->Cent); CHKERRQ(ierr); }
    if (user->lCent) { ierr = VecDestroy(&user->lCent); CHKERRQ(ierr); }
    if (user->GridSpace) { ierr = VecDestroy(&user->GridSpace); CHKERRQ(ierr); }
    if (user->lGridSpace) { ierr = VecDestroy(&user->lGridSpace); CHKERRQ(ierr); }
    if (user->Centx) { ierr = VecDestroy(&user->Centx); CHKERRQ(ierr); }
    if (user->Centy) { ierr = VecDestroy(&user->Centy); CHKERRQ(ierr); }
    if (user->Centz) { ierr = VecDestroy(&user->Centz); CHKERRQ(ierr); }
    if (user->lCentx) { ierr = VecDestroy(&user->lCentx); CHKERRQ(ierr); }
    if (user->lCenty) { ierr = VecDestroy(&user->lCenty); CHKERRQ(ierr); }
    if (user->lCentz) { ierr = VecDestroy(&user->lCentz); CHKERRQ(ierr); }
 
    // --- Group G: Turbulence Model Vectors (Finest level, conditional on les) ---
    if (user->Nu_t) { ierr = VecDestroy(&user->Nu_t); CHKERRQ(ierr); }
    if (user->lNu_t) { ierr = VecDestroy(&user->lNu_t); CHKERRQ(ierr); }
    if (user->CS) { ierr = VecDestroy(&user->CS); CHKERRQ(ierr); }
    if (user->lCs) { ierr = VecDestroy(&user->lCs); CHKERRQ(ierr); }
    if (user->Nu_Wall) { ierr = VecDestroy(&user->Nu_Wall); CHKERRQ(ierr); }
    if (user->lNu_Wall) { ierr = VecDestroy(&user->lNu_Wall); CHKERRQ(ierr); }
    if (user->Friction_Velocity) { ierr = VecDestroy(&user->Friction_Velocity); CHKERRQ(ierr); }
    if (user->lFriction_Velocity) { ierr = VecDestroy(&user->lFriction_Velocity); CHKERRQ(ierr); }
 
    // --- Group H: Particle Vectors (Finest level, conditional on np > 0) ---
    if (user->ParticleCount) { ierr = VecDestroy(&user->ParticleCount); CHKERRQ(ierr); }
    if (user->lParticleCount) { ierr = VecDestroy(&user->lParticleCount); CHKERRQ(ierr); }
    if (user->Psi) { ierr = VecDestroy(&user->Psi); CHKERRQ(ierr); }
    if (user->lPsi) { ierr = VecDestroy(&user->lPsi); CHKERRQ(ierr); }
 
    // --- Group I: Boundary Condition Vectors (All levels) ---
    if (user->Bcs.Ubcs) { ierr = VecDestroy(&user->Bcs.Ubcs); CHKERRQ(ierr); }
    if (user->Bcs.Uch) { ierr = VecDestroy(&user->Bcs.Uch); CHKERRQ(ierr); }
 
    // --- Group J: Post-Processing Vectors (Finest level, postprocessor mode) ---
    if (user->P_nodal) { ierr = VecDestroy(&user->P_nodal); CHKERRQ(ierr); }
    if (user->Ucat_nodal) { ierr = VecDestroy(&user->Ucat_nodal); CHKERRQ(ierr); }
    if (user->Qcrit) { ierr = VecDestroy(&user->Qcrit); CHKERRQ(ierr); }
    if (user->lQcrit) { ierr = VecDestroy(&user->lQcrit); CHKERRQ(ierr); }
    if (user->Qcrit_nodal) { ierr = VecDestroy(&user->Qcrit_nodal); CHKERRQ(ierr); }
    if (user->Psi_nodal) { ierr = VecDestroy(&user->Psi_nodal); CHKERRQ(ierr); }
 
    // --- Group K: Interpolation Vectors (Lazy allocation) ---
    if (user->fieldStatisticsStorage) {
        SimCtx *stats_ctx = user->simCtx;
        const PetscInt window_count = stats_ctx ? stats_ctx->fieldStatisticsWindowCount : 0;
        for (PetscInt w = 0; w < window_count; ++w) {
            ierr = PicurvWindowStorageDestroy(&user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
        }
        ierr = PetscFree(user->fieldStatisticsStorage); CHKERRQ(ierr);
        user->fieldStatisticsStorage = NULL;
    }
    if (user->CellScalarAtCorner) { ierr = VecDestroy(&user->CellScalarAtCorner); CHKERRQ(ierr); }
    if (user->lCellScalarAtCorner) { ierr = VecDestroy(&user->lCellScalarAtCorner); CHKERRQ(ierr); }
    if (user->CellVectorAtCorner) { ierr = VecDestroy(&user->CellVectorAtCorner); CHKERRQ(ierr); }
    if (user->PostScalar) { ierr = VecDestroy(&user->PostScalar); CHKERRQ(ierr); }
    if (user->lPostScalar) { ierr = VecDestroy(&user->lPostScalar); CHKERRQ(ierr); }
    if (user->PostScalarNodal) { ierr = VecDestroy(&user->PostScalarNodal); CHKERRQ(ierr); }
    if (user->PostVector) { ierr = VecDestroy(&user->PostVector); CHKERRQ(ierr); }
    if (user->lPostVector) { ierr = VecDestroy(&user->lPostVector); CHKERRQ(ierr); }
    if (user->PostVectorNodal) { ierr = VecDestroy(&user->PostVectorNodal); CHKERRQ(ierr); }
    if (user->lCellVectorAtCorner) { ierr = VecDestroy(&user->lCellVectorAtCorner); CHKERRQ(ierr); }
 
    // --- Group L: Implicit Solver Temporary Vectors (Destroyed after use, but check anyway) ---
    if (user->Rhs) { ierr = VecDestroy(&user->Rhs); CHKERRQ(ierr); }
    if (user->dUcont) { ierr = VecDestroy(&user->dUcont); CHKERRQ(ierr); }
    if (user->pUcont) { ierr = VecDestroy(&user->pUcont); CHKERRQ(ierr); }
 
    // --- Group M: Poisson Solver Vectors (kept for the run by PoissonSolver_Multigrid) ---
    if (user->B) { ierr = VecDestroy(&user->B); CHKERRQ(ierr); }
    if (user->R) { ierr = VecDestroy(&user->R); CHKERRQ(ierr); }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "All vectors destroyed for UserCtx.\n");
    PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "DestroyUserContext"
/**
 * @brief Internal helper implementation: `DestroyUserContext()`.
 * @details Local to this translation unit.
 */
PetscErrorCode DestroyUserContext(UserCtx *user)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!user) {
        LOG_ALLOW(LOCAL, LOG_WARNING, "DestroyUserContext called with NULL user pointer.\n");
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Destroying UserCtx at level %d...\n", user->thislevel);
 
    // --- Step 1: Destroy Boundary Condition System ---
    // This handles all BC handlers and their private data.
    ierr = BoundarySystem_Destroy(user); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "  Boundary system destroyed.\n");
 
    // --- Step 2: Destroy All Vectors ---
    // Handles ~74 Vec objects with proper NULL checking.
    ierr = DestroyUserVectors(user); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "  All vectors destroyed.\n");
 
    // --- Step 3: Destroy Matrix and Solver Objects ---
    // Destroy pressure-Poisson matrices and solver.
    if (user->A) {
        ierr = MatDestroy(&user->A); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix A destroyed.\n");
    }
    if (user->MR) {
        ierr = MatDestroy(&user->MR); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix MR destroyed.\n");
    }
    if (user->MP) {
        ierr = MatDestroy(&user->MP); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix MP destroyed.\n");
    }
    if (user->ksp) {
        ierr = KSPDestroy(&user->ksp); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  KSP solver destroyed.\n");
    }
    if (user->nullsp) {
        ierr = MatNullSpaceDestroy(&user->nullsp); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  MatNullSpace destroyed.\n");
    }
 
    // --- Step 4: Destroy Application Ordering ---
    if (user->ao) {
        ierr = AODestroy(&user->ao); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  AO destroyed.\n");
    }
 
    // --- Step 5: Destroy DM Objects ---
    // Destroy in reverse order of dependency: post_swarm, swarm, fda6, fda, da
    if (user->post_swarm) {
        ierr = DMDestroy(&user->post_swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  post_swarm DM destroyed.\n");
    }
    if (user->swarm) {
        ierr = DMDestroy(&user->swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  swarm DM destroyed.\n");
    }
    if (user->fda6) {
        ierr = DMDestroy(&user->fda6); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  fda6 DM destroyed.\n");
    }
    if (user->da) {
        ierr = DMDestroy(&user->da); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  da DM destroyed.\n");
    }
 
    // --- Step 6: Free PetscMalloc'd Arrays ---
    // Free arrays allocated with PetscMalloc1
    if (user->RankCellInfoMap) {
        ierr = PetscFree(user->RankCellInfoMap); CHKERRQ(ierr);
        user->RankCellInfoMap = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  RankCellInfoMap freed.\n");
    }
 
    LOG_ALLOW(LOCAL, LOG_INFO, "UserCtx at level %d fully destroyed.\n", user->thislevel);
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "FinalizeSimulation"
/**
 * @brief Implementation of \ref FinalizeSimulation().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/setup.h`.
 * @see FinalizeSimulation()
 */
PetscErrorCode FinalizeSimulation(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!simCtx) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "FinalizeSimulation called with NULL SimCtx pointer.\n");
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning simulation memory cleanup...\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
 
    // ============================================================================
    // PHASE 1: DESTROY MULTIGRID HIERARCHY (All UserCtx structures)
    // ============================================================================
 
    ierr = DestroySolutionConvergenceState(simCtx); CHKERRQ(ierr);
    ierr = DestroyFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
    ierr = ParticleFieldPlanDestroy(&simCtx->particleFieldPlan); CHKERRQ(ierr);
 
    if (simCtx->usermg.mgctx) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying multigrid hierarchy (%d levels)...\n",
                  simCtx->usermg.mglevels);
 
        // Destroy each UserCtx from finest to coarsest (reverse order is safer)
        for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
            UserCtx *user = simCtx->usermg.mgctx[level].user;
            if (user) {
                LOG_ALLOW(LOCAL, LOG_INFO, "  Destroying level %d of %d...\n",
                          level, simCtx->usermg.mglevels - 1);
                ierr = DestroyUserContext(user); CHKERRQ(ierr);
 
                // Free the UserCtx structure itself
                ierr = PetscFree(user); CHKERRQ(ierr);
                simCtx->usermg.mgctx[level].user = NULL;
            }
 
            // Destroy the MGCtx-level packer DM
            if (simCtx->usermg.mgctx[level].packer) {
                ierr = DMDestroy(&simCtx->usermg.mgctx[level].packer); CHKERRQ(ierr);
                LOG_ALLOW(LOCAL, LOG_DEBUG, "  MGCtx[%d].packer destroyed.\n", level);
            }
        }
 
        // Free the MGCtx array itself
        ierr = PetscFree(simCtx->usermg.mgctx); CHKERRQ(ierr);
        simCtx->usermg.mgctx = NULL;
        LOG_ALLOW(GLOBAL, LOG_INFO, "All multigrid levels destroyed.\n");
    }
 
    // ============================================================================
    // PHASE 2: DESTROY USERMG-LEVEL OBJECTS
    // ============================================================================
 
    if (simCtx->usermg.packer) {
        ierr = DMDestroy(&simCtx->usermg.packer); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.packer DM destroyed.\n");
    }
 
    if (simCtx->usermg.snespacker) {
        ierr = SNESDestroy(&simCtx->usermg.snespacker); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.snespacker SNES destroyed.\n");
    }
 
    // ============================================================================
    // PHASE 3: DESTROY SIMCTX-LEVEL OBJECTS
    // ============================================================================
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying SimCtx-level objects...\n");
 
    // --- PetscViewer for logging ---
    if (simCtx->logviewer) {
        ierr = PetscViewerDestroy(&simCtx->logviewer); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  logviewer destroyed.\n");
    }
 
    // --- Particle System DM ---
    if (simCtx->dm_swarm) {
        ierr = DMDestroy(&simCtx->dm_swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  dm_swarm destroyed.\n");
    }
 
    // --- BoundingBox List (Array of BoundingBox structs) ---
    if (simCtx->bboxlist) {
        ierr = PetscFree(simCtx->bboxlist); CHKERRQ(ierr);
        simCtx->bboxlist = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  bboxlist freed.\n");
    }
 
    // --- Boundary Condition Files (Array of strings) ---
    if (simCtx->bcs_files) {
        for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
            if (simCtx->bcs_files[i]) {
                ierr = PetscFree(simCtx->bcs_files[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
        simCtx->bcs_files = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  bcs_files array freed (%d files).\n", simCtx->num_bcs_files);
    }
 
    // --- Brownian Motion RNG ---
    if (simCtx->BrownianMotionRNG) {
        ierr = PetscRandomDestroy(&simCtx->BrownianMotionRNG); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  BrownianMotionRNG destroyed.\n");
    }
    // --- Post-Processing Parameters ---
    // pps is allocated with PetscNew and contains only static char arrays and basic types.
    // No internal dynamic allocations need to be freed.
    if (simCtx->pps) {
        ierr = PetscFree(simCtx->pps); CHKERRQ(ierr);
        simCtx->pps = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  PostProcessParams freed.\n");
    }
 
    // --- IBM/FSI Objects ---
    // Note: These are initialized to NULL and currently have no dedicated destroy functions.
    // If these modules are extended with cleanup routines, call them here.
    if (simCtx->ibm != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->ibm is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
    if (simCtx->ibmv != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->ibmv is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
    if (simCtx->fsi != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->fsi is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
 
    // --- Logging Allowed Functions (Array of strings) ---
    // Note: The logging system maintains its own copy via set_allowed_functions(),
    // so freeing simCtx->allowedFuncs will NOT affect LOG_ALLOW functionality.
    if (simCtx->allowedFuncs) {
        for (PetscInt i = 0; i < simCtx->nAllowed; i++) {
            if (simCtx->allowedFuncs[i]) {
                ierr = PetscFree(simCtx->allowedFuncs[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->allowedFuncs); CHKERRQ(ierr);
        simCtx->allowedFuncs = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  allowedFuncs array freed (%d functions).\n", simCtx->nAllowed);
    }
 
    // --- Profiling Critical Functions (Array of strings) ---
    if (simCtx->profilingSelectedFuncs) {
        for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; i++) {
            if (simCtx->profilingSelectedFuncs[i]) {
                ierr = PetscFree(simCtx->profilingSelectedFuncs[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
        simCtx->profilingSelectedFuncs = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  profilingSelectedFuncs array freed (%d functions).\n", simCtx->nProfilingSelectedFuncs);
    }
 
    // ============================================================================
    // PHASE 4: FINAL SUMMARY
    // ============================================================================
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "Simulation cleanup completed successfully.\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "All PETSc objects have been destroyed.\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
 
    ierr = PetscFree(simCtx); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}