logging.c

Go to the documentation of this file.

// logging.c
#include "logging.h"
#include "statistics_window.h"
#include "statistics_accumulator.h"
#include "AnalyticalSolutions.h"
#include "particle_field_catalog.h"
#include "verification_sources.h"
 
/* Maximum temporary buffer size for converting numbers to strings */
#define TMP_BUF_SIZE 128
 
// --------------------- Static Variable for Log Level ---------------------
 
/**
 * @brief Static variable to cache the current logging level.
 *
 * Initialized to -1 to indicate that the log level has not been set yet.
 */
static LogLevel current_log_level = -1;
 
// --------------------- Static Variables for Allow-List -------------------
 
/**
 * @brief Global/static array of function names allowed to log.
 */
static char** gAllowedFunctions = NULL;
 
/**
 * @brief Number of entries in the gAllowedFunctions array.
 */
static int gNumAllowed = 0;
 
enum {
    SEARCH_METRIC_SUM_SEARCH_ATTEMPTS = 0,
    SEARCH_METRIC_SUM_SEARCH_POPULATION,
    SEARCH_METRIC_SUM_SEARCH_LOCATED,
    SEARCH_METRIC_SUM_SEARCH_LOST,
    SEARCH_METRIC_SUM_TRAVERSAL_STEPS,
    SEARCH_METRIC_SUM_RESEARCH,
    SEARCH_METRIC_SUM_TIE_BREAKS,
    SEARCH_METRIC_SUM_BOUNDARY_CLAMPS,
    SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS,
    SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK,
    SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS,
    SEARCH_METRIC_MAX_TRAVERSAL_STEPS,
    SEARCH_METRIC_MAX_PASS_DEPTH,
    SEARCH_METRIC_REDUCTION_LEN
};
 
/**
 * @brief Internal reduction callback for packed search metrics.
 * @details Local to this translation unit.
 */
static void SearchMetricsReduceOp(void *invec, void *inoutvec, int *len, MPI_Datatype *datatype)
{
    PetscReal *in = (PetscReal *)invec;
    PetscReal *inout = (PetscReal *)inoutvec;
    (void)datatype;
 
    for (int idx = 0; idx < *len; idx += SEARCH_METRIC_REDUCTION_LEN) {
        inout[idx + SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] += in[idx + SEARCH_METRIC_SUM_SEARCH_ATTEMPTS];
        inout[idx + SEARCH_METRIC_SUM_SEARCH_POPULATION] += in[idx + SEARCH_METRIC_SUM_SEARCH_POPULATION];
        inout[idx + SEARCH_METRIC_SUM_SEARCH_LOCATED] += in[idx + SEARCH_METRIC_SUM_SEARCH_LOCATED];
        inout[idx + SEARCH_METRIC_SUM_SEARCH_LOST] += in[idx + SEARCH_METRIC_SUM_SEARCH_LOST];
        inout[idx + SEARCH_METRIC_SUM_TRAVERSAL_STEPS] += in[idx + SEARCH_METRIC_SUM_TRAVERSAL_STEPS];
        inout[idx + SEARCH_METRIC_SUM_RESEARCH] += in[idx + SEARCH_METRIC_SUM_RESEARCH];
        inout[idx + SEARCH_METRIC_SUM_TIE_BREAKS] += in[idx + SEARCH_METRIC_SUM_TIE_BREAKS];
        inout[idx + SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] += in[idx + SEARCH_METRIC_SUM_BOUNDARY_CLAMPS];
        inout[idx + SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] += in[idx + SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS];
        inout[idx + SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] += in[idx + SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK];
        inout[idx + SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] += in[idx + SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS];
        inout[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS] = PetscMax(inout[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS],
                                                                   in[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS]);
        inout[idx + SEARCH_METRIC_MAX_PASS_DEPTH] = PetscMax(inout[idx + SEARCH_METRIC_MAX_PASS_DEPTH],
                                                             in[idx + SEARCH_METRIC_MAX_PASS_DEPTH]);
    }
}
 
// --------------------- Function Implementations ---------------------
 
/**
 * @brief Implementation of \ref get_log_level().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see get_log_level()
 */
LogLevel get_log_level() {
    if (current_log_level == -1) { // Log level not set yet
        const char *env = getenv("LOG_LEVEL");
        if (!env) {
            current_log_level = LOG_ERROR; // Default level
        }
        else if (strcmp(env, "DEBUG") == 0) {
            current_log_level = LOG_DEBUG;
        }
        else if (strcmp(env, "INFO") == 0) {
            current_log_level = LOG_INFO;
        }
        else if (strcmp(env, "WARNING") == 0) {
            current_log_level = LOG_WARNING;
        }
        else if (strcmp(env, "VERBOSE") == 0) {
            current_log_level = LOG_VERBOSE;
        }
        else if (strcmp(env, "TRACE") == 0) {
            current_log_level = LOG_TRACE;
        }
        else {
            current_log_level = LOG_ERROR; // Default if unrecognized
        }
    }
    return current_log_level;
}
 
/**
 * @brief Internal helper implementation: `print_log_level()`.
 * @details Local to this translation unit.
 */
PetscErrorCode print_log_level(void)
{
  PetscMPIInt     rank;
  PetscErrorCode  ierr;
  int             level;
  const char     *level_name;
 
  PetscFunctionBeginUser;
  /* get MPI rank */
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
 
  /* decide level name */
  level       = get_log_level();
  level_name = (level == LOG_ERROR)   ? "ERROR"   :
               (level == LOG_WARNING) ? "WARNING" :
               (level == LOG_INFO)    ? "INFO"    :
               (level == LOG_DEBUG)   ? "DEBUG"   :
               (level == LOG_VERBOSE) ? "VERBOSE" :
               (level == LOG_TRACE)   ? "TRACE"   :
               "UNKNOWN";
 
  /* print it out */
  ierr = PetscPrintf(PETSC_COMM_SELF,
                     "Current log level: %s (%d) | rank: %d\n",
                     level_name, level, (int)rank);
  CHKERRMPI(ierr);
 
  PetscFunctionReturn(PETSC_SUCCESS);
}
 
/**
 * @brief Implementation of \ref set_allowed_functions().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see set_allowed_functions()
 */
void set_allowed_functions(const char** functionList, int count)
{
    // 1. Free any existing entries
    if (gAllowedFunctions) {
        for (int i = 0; i < gNumAllowed; ++i) {
            free(gAllowedFunctions[i]); // each was strdup'ed
        }
        free(gAllowedFunctions);
        gAllowedFunctions = NULL;
        gNumAllowed = 0;
    }
 
    // 2. Allocate new array
    if (count > 0) {
        gAllowedFunctions = (char**)malloc(sizeof(char*) * count);
    }
 
    // 3. Copy the new entries
    for (int i = 0; i < count; ++i) {
        // strdup is a POSIX function. If not available, implement your own string copy.
        gAllowedFunctions[i] = strdup(functionList[i]);
    }
    gNumAllowed = count;
}
 
/**
 * @brief Implementation of \ref is_function_allowed().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see is_function_allowed()
 */
PetscBool is_function_allowed(const char* functionName)
{
    /* no list ⇒ allow all */
    if (gNumAllowed == 0) {
        return PETSC_TRUE;
    }
 
    /* otherwise only the listed functions are allowed */
    for (int i = 0; i < gNumAllowed; ++i) {
        if (strcmp(gAllowedFunctions[i], functionName) == 0) {
            return PETSC_TRUE;
        }
    }
    return PETSC_FALSE;
}
 
/**
 * @brief Implementation of \ref LOG_CELL_VERTICES().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_CELL_VERTICES()
 */
PetscErrorCode LOG_CELL_VERTICES(const Cell *cell, PetscMPIInt rank)
{
 
    // Validate input pointers
    if (cell == NULL) {
      LOG_ALLOW(LOCAL,LOG_ERROR, "'cell' is NULL.\n");
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "LOG_CELL_VERTICES - Input parameter 'cell' is NULL.");
    }
 
      LOG_ALLOW(LOCAL,LOG_VERBOSE, "Rank %d, Cell Vertices:\n", rank);
        for(int i = 0; i < 8; i++){ 
      LOG_ALLOW(LOCAL,LOG_VERBOSE, "  Vertex[%d]: (%.2f, %.2f, %.2f)\n", 
                       i, cell->vertices[i].x, cell->vertices[i].y, cell->vertices[i].z);
        }
 
    return 0; // Indicate successful execution
}
 
 
/**
 * @brief Implementation of \ref LOG_FACE_DISTANCES().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_FACE_DISTANCES()
 */
PetscErrorCode LOG_FACE_DISTANCES(PetscReal* d)
{
 
    // Validate input array
    if (d == NULL) {
      LOG_ALLOW(LOCAL,LOG_ERROR, "  'd' is NULL.\n");
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "  Input array 'd' is NULL.");
    }
 
    PetscPrintf(PETSC_COMM_SELF, "  Face Distances:\n");
    PetscPrintf(PETSC_COMM_SELF, "  LEFT(%d):   %.15f\n", LEFT, d[LEFT]);
    PetscPrintf(PETSC_COMM_SELF, "  RIGHT(%d):  %.15f\n", RIGHT, d[RIGHT]);
    PetscPrintf(PETSC_COMM_SELF, "  BOTTOM(%d): %.15f\n", BOTTOM, d[BOTTOM]);
    PetscPrintf(PETSC_COMM_SELF, "  TOP(%d):    %.15f\n", TOP, d[TOP]);
    PetscPrintf(PETSC_COMM_SELF, "  FRONT(%d):  %.15f\n", FRONT, d[FRONT]);
    PetscPrintf(PETSC_COMM_SELF, "  BACK(%d):   %.15f\n", BACK, d[BACK]);
 
    return 0; // Indicate successful execution
}
 
/*
 * Helper function: Converts an integer (of type int) to a string.
 */
static void IntToStr(int value, char *buf, size_t bufsize)
{
    snprintf(buf, bufsize, "%d", value);
}
 
/*
 * Helper function: Converts a 64‐bit integer to a string.
 */
static void Int64ToStr(PetscInt64 value, char *buf, size_t bufsize)
{
    snprintf(buf, bufsize, "%ld", value);
}
 
/*
 * Helper function: Converts three integers into a formatted string "(i, j, k)".
 */
static void CellToStr(const PetscInt *cell, char *buf, size_t bufsize)
{
    snprintf(buf, bufsize, "(%d, %d, %d)", cell[0], cell[1], cell[2]);
}
 
/*
 * Helper function: Converts three PetscReal values into a formatted string "(x, y, z)".
 */
static void TripleRealToStr(const PetscReal *arr, char *buf, size_t bufsize)
{
    snprintf(buf, bufsize, "(%.4f, %.4f, %.4f)", arr[0], arr[1], arr[2]);
}
 
/*
 * Helper function: Computes the maximum string length for each column (across all particles).
 *
 * The function examines every particle (from 0 to nParticles-1) and converts the value to a
 * string using the helper functions above. The maximum length is stored in the pointers provided.
 *
 * @param nParticles       Number of particles.
 * @param ranks            Array of particle MPI ranks.
 * @param pids             Array of particle IDs.
 * @param cellIDs          Array of cell IDs (stored consecutively, 3 per particle).
 * @param positions        Array of positions (3 per particle).
 * @param velocities       Array of velocities (3 per particle).
 * @param weights          Array of weights (3 per particle).
 * @param wRank            [out] Maximum width for Rank column.
 * @param wPID             [out] Maximum width for PID column.
 * @param wCell            [out] Maximum width for Cell column.
 * @param wPos             [out] Maximum width for Position column.
 * @param wVel             [out] Maximum width for Velocity column.
 * @param wWt              [out] Maximum width for Weights column.
 */
static PetscErrorCode ComputeMaxColumnWidths(PetscInt nParticles,
                                               const PetscMPIInt *ranks,
                                               const PetscInt64  *pids,
                                               const PetscInt  *cellIDs,
                                               const PetscReal   *positions,
                                               const PetscReal   *velocities,
                                               const PetscReal   *weights,
                                               int *wRank, int *wPID, int *wCell,
                                               int *wPos, int *wVel, int *wWt)
{
    char tmp[TMP_BUF_SIZE];
 
    *wRank = strlen("Rank");  /* Start with the header label lengths */
    *wPID  = strlen("PID");
    *wCell = strlen("Cell (i,j,k)");
    *wPos  = strlen("Position (x,y,z)");
    *wVel  = strlen("Velocity (x,y,z)");
    *wWt   = strlen("Weights (a1,a2,a3)");
 
    for (PetscInt i = 0; i < nParticles; i++) {
        /* Rank */
        IntToStr(ranks[i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wRank) *wRank = (int)strlen(tmp);
 
        /* PID */
        Int64ToStr(pids[i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wPID) *wPID = (int)strlen(tmp);
 
        /* Cell: use the three consecutive values */
        CellToStr(&cellIDs[3 * i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wCell) *wCell = (int)strlen(tmp);
 
        /* Position */
        TripleRealToStr(&positions[3 * i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wPos) *wPos = (int)strlen(tmp);
 
        /* Velocity */
        TripleRealToStr(&velocities[3 * i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wVel) *wVel = (int)strlen(tmp);
 
        /* Weights */
        TripleRealToStr(&weights[3 * i], tmp, TMP_BUF_SIZE);
        if ((int)strlen(tmp) > *wWt) *wWt = (int)strlen(tmp);
    }
    return 0;
}
 
/*
 * Helper function: Builds a format string for a table row.
 *
 * The format string will include proper width specifiers for each column.
 * For example, it might create something like:
 *
 * "| %-6s | %-8s | %-20s | %-25s | %-25s | %-25s |\n"
 *
 * @param wRank   Maximum width for the Rank column.
 * @param wPID    Maximum width for the PID column.
 * @param wCell   Maximum width for the Cell column.
 * @param wPos    Maximum width for the Position column.
 * @param wVel    Maximum width for the Velocity column.
 * @param wWt     Maximum width for the Weights column.
 * @param fmtStr  Buffer in which to build the format string.
 * @param bufSize Size of fmtStr.
 */
static void BuildRowFormatString(PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt, char *fmtStr, size_t bufSize)
{
    // Build a format string using snprintf.
    // We assume that the Rank is an int (%d), PID is a 64-bit int (%ld)
    // and the remaining columns are strings (which have been formatted already).
    snprintf(fmtStr, bufSize,
             "| %%-%dd | %%-%dd | %%-%ds | %%-%ds | %%-%ds | %%-%ds |\n",
             wRank, wPID, wCell, wPos, wVel, wWt);
}
 
/*
 * Helper function: Builds a header string for the table using column titles.
 */
static void BuildHeaderString(char *headerStr, size_t bufSize, PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt)
{
    snprintf(headerStr, bufSize,
             "| %-*s | %-*s | %-*s | %-*s | %-*s | %-*s |\n",
             (int)wRank, "Rank",
             (int)wPID, "PID",
             (int)wCell, "Cell (i,j,k)",
             (int)wPos, "Position (x,y,z)",
             (int)wVel, "Velocity (x,y,z)",
             (int)wWt, "Weights (a1,a2,a3)");
}
 
/**
 * @brief Implementation of \ref LOG_PARTICLE_FIELDS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_PARTICLE_FIELDS()
 */
PetscErrorCode LOG_PARTICLE_FIELDS(UserCtx* user, PetscInt printInterval)
{
    DM swarm = user->swarm;
    PetscErrorCode ierr;
    PetscInt localNumParticles;
    PetscReal *positions = NULL;
    PetscInt64 *particleIDs = NULL;
    PetscMPIInt *particleRanks = NULL;
    PetscInt *cellIDs = NULL;
    PetscReal *weights = NULL;
    PetscReal *velocities = NULL;
    PetscMPIInt rank;
    
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_INFO, "Rank %d is retrieving particle data.\n", rank);
 
    ierr = DMSwarmGetLocalSize(swarm, &localNumParticles); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_DEBUG,"Rank %d has %d particles.\n", rank, localNumParticles);
 
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&positions); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void**)&particleIDs); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_RANK), NULL, NULL, (void**)&particleRanks); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT), NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), NULL, NULL, (void**)&velocities); CHKERRQ(ierr);
    
    /* Compute maximum column widths. */
    int wRank, wPID, wCell, wPos, wVel, wWt;
    wRank = wPID = wCell = wPos = wVel = wWt = 0;
    ierr = ComputeMaxColumnWidths(localNumParticles, particleRanks, particleIDs, cellIDs,
                                  positions, velocities, weights,
                                  &wRank, &wPID, &wCell, &wPos, &wVel, &wWt); CHKERRQ(ierr);
 
    /* Build a header string and a row format string. */
    char headerFmt[256];
    char rowFmt[256];
    BuildHeaderString(headerFmt, sizeof(headerFmt), wRank, wPID, wCell, wPos, wVel, wWt);
    BuildRowFormatString(wRank, wPID, wCell, wPos, wVel, wWt, rowFmt, sizeof(rowFmt));
    
    /* Print header (using synchronized printing for parallel output). */
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", headerFmt); CHKERRQ(ierr);
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
    
    /* Loop over particles and print every printInterval-th row. */
    char rowStr[256];
    for (PetscInt i = 0; i < localNumParticles; i++) {
        if (i % printInterval == 0) {
      // ------- DEBUG 
      //char cellStr[TMP_BUF_SIZE], posStr[TMP_BUF_SIZE], velStr[TMP_BUF_SIZE], wtStr[TMP_BUF_SIZE];
      //CellToStr(&cellIDs[3*i], cellStr, TMP_BUF_SIZE);
      //TripleRealToStr(&positions[3*i], posStr, TMP_BUF_SIZE);
      //TripleRealToStr(&velocities[3*i], velStr, TMP_BUF_SIZE);
      // TripleRealToStr(&weights[3*i], wtStr, TMP_BUF_SIZE);
            
      // if (rank == 0) { // Or whatever rank is Rank 0
        //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Particle %lld: PID=%lld, Rank=%d\n", (long long)i, (long long)particleIDs[i], particleRanks[i]);
        //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Raw Pos: (%.10e, %.10e, %.10e)\n", positions[3*i+0], positions[3*i+1], positions[3*i+2]);
        //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Str Pos: %s\n", posStr);
        //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Raw Vel: (%.10e, %.10e, %.10e)\n", velocities[3*i+0], velocities[3*i+1], velocities[3*i+2]);
        // PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Str Vel: %s\n", velStr);
        // Add similar for cell, weights
        // PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] About to build rowStr for particle %lld\n", (long long)i);
        //  fflush(stdout);
        // }
 
      //  snprintf(rowStr, sizeof(rowStr), rowFmt,
          //           particleRanks[i],
          //           particleIDs[i],
          //           cellStr,
          //           posStr,
          //           velStr,
      //         wtStr);
 
         
         //    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr);
      
         // ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr); 
      
      // -------- DEBUG
            /* Format the row by converting each field to a string first.
             * We use temporary buffers and then build the row string.
             */
      
       char cellStr[TMP_BUF_SIZE], posStr[TMP_BUF_SIZE], velStr[TMP_BUF_SIZE], wtStr[TMP_BUF_SIZE];
            CellToStr(&cellIDs[3*i], cellStr, TMP_BUF_SIZE);
            TripleRealToStr(&positions[3*i], posStr, TMP_BUF_SIZE);
            TripleRealToStr(&velocities[3*i], velStr, TMP_BUF_SIZE);
            TripleRealToStr(&weights[3*i], wtStr, TMP_BUF_SIZE);
            
            /* Build the row string. Note that for the integer fields we can use the row format string. */
            snprintf(rowStr, sizeof(rowStr), rowFmt,
                     particleRanks[i],
                     particleIDs[i],
                     cellStr,
                     posStr,
                     velStr,
                     wtStr);
     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr);
        }
    }
 
 
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"); CHKERRQ(ierr);
    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_DEBUG,"Completed printing on Rank %d.\n", rank);
 
    /* Restore fields */
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&positions); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void**)&particleIDs); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_RANK), NULL, NULL, (void**)&particleRanks); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT), NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), NULL, NULL, (void**)&velocities); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL,LOG_DEBUG, "Restored all particle fields.\n");
    return 0;
}
 
/**
 * @brief Implementation of \ref IsParticleConsoleSnapshotEnabled().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see IsParticleConsoleSnapshotEnabled()
 */
 
PetscBool IsParticleConsoleSnapshotEnabled(const SimCtx *simCtx)
{
    if (!simCtx) {
        return PETSC_FALSE;
    }
    return (PetscBool)(simCtx->np > 0 &&
                       simCtx->particleConsoleOutputFreq > 0 &&
                       get_log_level() >= LOG_INFO);
}
 
/**
 * @brief Implementation of \ref ShouldEmitPeriodicParticleConsoleSnapshot().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ShouldEmitPeriodicParticleConsoleSnapshot()
 */
 
PetscBool ShouldEmitPeriodicParticleConsoleSnapshot(const SimCtx *simCtx, PetscInt completed_step)
{
    return (PetscBool)(IsParticleConsoleSnapshotEnabled(simCtx) &&
                       completed_step > 0 &&
                       completed_step % simCtx->particleConsoleOutputFreq == 0);
}
 
/**
 * @brief Implementation of \ref EmitParticleConsoleSnapshot().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see EmitParticleConsoleSnapshot()
 */
 
PetscErrorCode EmitParticleConsoleSnapshot(UserCtx *user, SimCtx *simCtx, PetscInt step)
{
    PetscErrorCode ierr;
 
    PetscFunctionBeginUser;
    LOG(GLOBAL, LOG_INFO, "Particle states at step %d:\n", step);
    ierr = LOG_PARTICLE_FIELDS(user, simCtx->LoggingFrequency); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}
 
 
/**
 * @brief Remove leading and trailing whitespace from a mutable configuration string.
 */
static void trim(char *s)
{
    if (!s) return;
 
    /* ---- 1. strip leading blanks ----------------------------------- */
    char *p = s;
    while (*p && isspace((unsigned char)*p))
        ++p;
 
    if (p != s)                      /* move the trimmed text forward   */
        memmove(s, p, strlen(p) + 1);   /* +1 to copy the final NUL     */
 
    /* ---- 2. strip trailing blanks ---------------------------------- */
    size_t len = strlen(s);
    while (len > 0 && isspace((unsigned char)s[len - 1]))
        s[--len] = '\0';
}
 
/* ------------------------------------------------------------------------- */
/**
 * @brief Implementation of \ref LoadAllowedFunctionsFromFile().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LoadAllowedFunctionsFromFile()
 */
PetscErrorCode LoadAllowedFunctionsFromFile(const char   filename[],
                                            char      ***funcsOut,
                                            PetscInt   *nOut)
{
  FILE          *fp    = NULL;
  char         **funcs = NULL;
  size_t         cap   = 16;   /* initial capacity */
  size_t         n     = 0;    /* number of names  */
  char           line[PETSC_MAX_PATH_LEN];
  PetscErrorCode ierr;
 
  PetscFunctionBegin;
 
  /* ---------------------------------------------------------------------- */
  /* 1. Open file                                                           */
  fp = fopen(filename, "r");
  if (!fp) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
                    "Cannot open %s", filename);
 
  /* 2. Allocate initial pointer array                                      */
  ierr = PetscMalloc1(cap, &funcs); CHKERRQ(ierr);
 
  /* 3. Read file line by line                                              */
  while (fgets(line, sizeof line, fp)) {
    /* Strip everything after a comment character '#'. */
    char *hash = strchr(line, '#');
    if (hash) *hash = '\0';
 
    trim(line);                 /* remove leading/trailing blanks */
    if (!*line) continue;       /* skip if empty                  */
 
    /* Grow the array if necessary */
    if (n == cap) {
      cap *= 2;
      ierr = PetscRealloc(cap * sizeof(*funcs), (void **)&funcs); CHKERRQ(ierr);
    }
 
    /* Deep‑copy the cleaned identifier */
    ierr = PetscStrallocpy(line, &funcs[n++]); CHKERRQ(ierr);
  }
  fclose(fp);
 
  /* 4. Return results to caller                                           */
  *funcsOut = funcs;
  *nOut     = (PetscInt)n;
 
  PetscFunctionReturn(0);
}
 
/* ------------------------------------------------------------------------- */
/**
 * @brief Internal helper implementation: `FreeAllowedFunctions()`.
 * @details Local to this translation unit.
 */
PetscErrorCode FreeAllowedFunctions(char **funcs, PetscInt n)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  if (funcs) {
    for (PetscInt i = 0; i < n; ++i) {
      ierr = PetscFree(funcs[i]); CHKERRQ(ierr);
    }
    ierr = PetscFree(funcs); CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref BCFaceToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see BCFaceToString()
 */
const char* BCFaceToString(BCFace face) {
    switch (face) {
        case BC_FACE_NEG_X: return "-Xi (I-Min)";
        case BC_FACE_POS_X: return "+Xi (I-Max)";
        case BC_FACE_NEG_Y: return "-Eta (J-Min)";
        case BC_FACE_POS_Y: return "+Eta (J-Max)";
        case BC_FACE_NEG_Z: return "-Zeta (K-Min)";
        case BC_FACE_POS_Z: return "+Zeta (K-Max)";
        default:            return "Unknown Face";
    }
}
 
/**
 * @brief Implementation of \ref InitialConditionModeToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see InitialConditionModeToString()
 */
const char* InitialConditionModeToString(InitialConditionMode mode)
{
    switch(mode){
        case IC_MODE_ZERO: return "Zero";
        case IC_MODE_CONSTANT_CARTESIAN: return "Cartesian Constant";
        case IC_MODE_POISEUILLE: return "Poiseuille";
        case IC_MODE_CONSTANT_STREAMWISE: return "Streamwise Constant";
        case IC_MODE_FILE: return "File";
        default: return "Unknown Initial Condition";
    }
}
 
/*
 * Converts a FlowDirection enum value to its token string. The public header
 * owns the rendered API contract.
 */
const char* FlowDirectionToString(FlowDirection fd)
{
    switch ((int)fd) {
        case FLOW_DIR_POS_XI:   return "+Xi";
        case FLOW_DIR_NEG_XI:   return "-Xi";
        case FLOW_DIR_POS_ETA:  return "+Eta";
        case FLOW_DIR_NEG_ETA:  return "-Eta";
        case FLOW_DIR_POS_ZETA: return "+Zeta";
        case FLOW_DIR_NEG_ZETA: return "-Zeta";
        default:                return "from INLET";
    }
}
 
/**
 * @brief Implementation of \ref ParticleInitializationToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ParticleInitializationToString()
 */
const char* ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
{
    switch(ParticleInitialization){
        case PARTICLE_INIT_SURFACE_RANDOM: return "Surface: Random";
        case PARTICLE_INIT_VOLUME:         return "Volume";
        case PARTICLE_INIT_POINT_SOURCE:   return "Point Source";
        case PARTICLE_INIT_SURFACE_EDGES:  return "Surface: At edges";
        default: return "Unknown Particle Initialization";
    }
}
 
/**
 * @brief Implementation of \ref LESModelToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LESModelToString()
 */
const char* LESModelToString(LESModelType LESFlag)
{
    switch(LESFlag){
        case NO_LES_MODEL: return "No LES";
        case CONSTANT_SMAGORINSKY: return "Constant Smagorinsky";
        case DYNAMIC_SMAGORINSKY: return "Dynamic Smagorinsky";
        case VREMAN: return "Vreman";
        case WALE: return "WALE";
        default: return "Unknown LES Flag";
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "LESFilterWidthModelToString"
/**
 * @brief Implementation of \ref LESFilterWidthModelToString().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 * @see LESFilterWidthModelToString()
 */
const char* LESFilterWidthModelToString(LESFilterWidthModel model)
{
    /* These spellings are the user-facing YAML values, so a banner line can be pasted
       back into a case file without translation. */
    switch(model){
        case LES_FILTER_WIDTH_CUBE_ROOT_VOLUME: return "cube_root_volume";
        case LES_FILTER_WIDTH_GEOMETRIC_MEAN:   return "geometric_mean";
        case LES_FILTER_WIDTH_MAX_EDGE:         return "max_edge";
        case LES_FILTER_WIDTH_SCOTTI:           return "scotti";
        default: return "unknown";
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "LESTestFilterKernelToString"
/**
 * @brief Implementation of \ref LESTestFilterKernelToString().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 * @see LESTestFilterKernelToString()
 */
const char* LESTestFilterKernelToString(LESTestFilterKernel kernel)
{
    switch(kernel){
        case LES_TEST_FILTER_VOLUME_WEIGHTED_BOX: return "volume_weighted_box";
        case LES_TEST_FILTER_SIMPSON_IK:          return "simpson_ik";
        default: return "unknown";
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "LESAveragingModeToString"
/**
 * @brief Implementation of \ref LESAveragingModeToString().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 * @see LESAveragingModeToString()
 */
const char* LESAveragingModeToString(LESAveragingMode mode)
{
    switch(mode){
        case LES_AVERAGING_LOCAL:       return "local";
        case LES_AVERAGING_HOMOGENEOUS: return "homogeneous";
        case LES_AVERAGING_GLOBAL:      return "global";
        default: return "unknown";
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "LESClipModeToString"
/**
 * @brief Implementation of \ref LESClipModeToString().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 * @see LESClipModeToString()
 */
const char* LESClipModeToString(LESClipMode mode)
{
    switch(mode){
        case LES_CLIP_CLAMP:         return "clamp";
        case LES_CLIP_CLIP_NEGATIVE: return "clip_negative";
        case LES_CLIP_NONE:          return "none";
        default: return "unknown";
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "WallFunctionModelToString"
/**
 * @brief Implementation of \ref WallFunctionModelToString().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 * @see WallFunctionModelToString()
 */
const char* WallFunctionModelToString(WallFunctionModel model)
{
    /* User-facing YAML spellings, so a banner line can be pasted back into a case. */
    switch(model){
        case WALL_FUNCTION_NONE:    return "none";
        case WALL_FUNCTION_LOG_LAW: return "log_law";
        case WALL_FUNCTION_WERNER:  return "werner";
        case WALL_FUNCTION_CABOT:   return "cabot";
        default: return "unknown";
    }
}
 
/**
 * @brief Implementation of \ref MomentumSolverTypeToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see MomentumSolverTypeToString()
 */
const char* MomentumSolverTypeToString(MomentumSolverType SolverFlag)
{
    switch(SolverFlag){
        case MOMENTUM_SOLVER_EXPLICIT_RK: return "Explicit 4 stage Runge-Kutta ";
        case MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK: return "Dual Time Picard with 4-stage Jameson RK Smoothing";
        case MOMENTUM_SOLVER_NEWTON_KRYLOV: return "Newton Krylov";
        default: return "Unknown Momentum Solver Type";
    }
}
 
/**
 * @brief Implementation of \ref BCTypeToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see BCTypeToString()
 */
const char* BCTypeToString(BCType type) {
    switch (type) {
      //  case DIRICHLET: return "DIRICHLET";
      //  case NEUMANN:   return "NEUMANN";
        case WALL:      return "WALL";
        case INLET:     return "INLET";
        case OUTLET:    return "OUTLET";
        case FARFIELD:  return "FARFIELD";
        case PERIODIC:  return "PERIODIC";
        case INTERFACE: return "INTERFACE";
 
    // case CUSTOM:    return "CUSTOM";
        default:        return "Unknown BC Type";
    }
}
 
/**
 * @brief Internal helper implementation: `BCHandlerTypeToString()`.
 * @details Local to this translation unit.
 */
const char* BCHandlerTypeToString(BCHandlerType handler_type) {
    switch (handler_type) {
        // Wall & Symmetry Handlers
        case BC_HANDLER_WALL_NOSLIP:             return "noslip";
        case BC_HANDLER_WALL_MOVING:             return "moving";
        case BC_HANDLER_SYMMETRY_PLANE:          return "symmetry_plane";
 
        // Inlet Handlers
        case BC_HANDLER_INLET_CONSTANT_VELOCITY: return "constant_velocity";
        case BC_HANDLER_INLET_PARABOLIC:         return "parabolic";
        case BC_HANDLER_INLET_PROFILE_FROM_FILE: return "prescribed_flow";
 
        // Outlet Handlers
        case BC_HANDLER_OUTLET_CONSERVATION:     return "conservation";
        case BC_HANDLER_OUTLET_PRESSURE:         return "pressure";
 
        // Other Physical Handlers
        case BC_HANDLER_FARFIELD_NONREFLECTING:  return "nonreflecting";
 
        // Multi-Block / Interface Handlers
        case BC_HANDLER_PERIODIC_GEOMETRIC:      return "geometric";
        case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX:   return "constant flux";
        case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX:    return "initial flux";
        case BC_HANDLER_INTERFACE_OVERSET:       return "overset";
 
        // Default case
        case BC_HANDLER_UNDEFINED:
        default:                                 return "UNKNOWN_HANDLER";
    }
}
 
/**
 * @brief Implementation of \ref DualMonitorDestroy().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see DualMonitorDestroy()
 */
PetscErrorCode DualMonitorDestroy(void **ctx)
{
    DualMonitorCtx *monctx = (DualMonitorCtx*)*ctx;
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
    if(!rank && monctx->file_handle){
      fclose(monctx->file_handle);
    }
    
    ierr = PetscFree(monctx); CHKERRQ(ierr);
    *ctx = NULL;
    PetscFunctionReturn(0);
}
 
/**
 * @brief A custom KSP monitor that logs the true residual to a file and optionally to the console.
 *
 * This function replicates the behavior of KSPMonitorTrueResidualNorm by calculating
 * the true residual norm ||b - Ax|| itself. It unconditionally logs to a file
 * viewer and conditionally logs to the console based on a flag in the context.
 *
 * @param ksp   The Krylov subspace context.
 * @param it    The current iteration number.
 * @param rnorm The preconditioned residual norm (ignored, we compute our own).
 * @param ctx   A pointer to the DualMonitorCtx structure.
 * @return      PetscErrorCode 0 on success.
 */
#undef __FUNCT__
#define __FUNCT__ "DualKSPMonitor"
/**
 * @brief Implementation of \ref DualKSPMonitor().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see DualKSPMonitor()
 */
 
PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
{
    DualMonitorCtx *monctx = (DualMonitorCtx*)ctx;
    PetscErrorCode ierr;
    PetscReal      trnorm, relnorm;
    Vec            r;
    char           norm_buf[256];
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
 
    // 1. Calculate the true residual norm.
    ierr = KSPBuildResidual(ksp, NULL, NULL, &r); CHKERRQ(ierr);
    ierr = VecNorm(r, NORM_2, &trnorm); CHKERRQ(ierr);
    ierr = VecDestroy(&r); CHKERRQ(ierr);
 
    // 2. On the first iteration, compute and store the norm of the RHS vector `b`.
    if (it == 0) {
        Vec b;
        ierr = KSPGetRhs(ksp, &b); CHKERRQ(ierr);
        ierr = VecNorm(b, NORM_2, &monctx->bnorm); CHKERRQ(ierr);
    }
 
    if(!rank){
    // 3. Compute the relative norm and format the output string.
    if (monctx->bnorm > 1.e-15) {
      relnorm = trnorm / monctx->bnorm;
      sprintf(norm_buf, "ts: %-5d | block: %-2d | iter: %-3d | Unprecond Norm: %12.5e | True Norm: %12.5e | Rel Norm: %12.5e",(int)monctx->step, (int)monctx->block_id, (int)it, (double)rnorm, (double)trnorm, (double)relnorm);
    } else {
      sprintf(norm_buf,"ts: %-5d | block: %-2d | iter: %-3d | Unprecond Norm: %12.5e | True Norm: %12.5e",(int)monctx->step, (int)monctx->block_id, (int)it, (double)rnorm, (double)trnorm);
    }
 
    // 4. Log to the file viewer (unconditionally).
    if(monctx->file_handle){
      ierr = PetscFPrintf(PETSC_COMM_SELF,monctx->file_handle,"%s\n", norm_buf); CHKERRQ(ierr);
    }
    // 5. Log to the console (conditionally).
    if (monctx->log_to_console) {
      PetscFPrintf(PETSC_COMM_SELF,stdout, "%s\n", norm_buf); CHKERRQ(ierr);
    }
 
    } //rank
 
    PetscFunctionReturn(0);
}
 
#define SOLUTION_CONVERGENCE_FLUID_THRESHOLD 0.1
#define SOLUTION_CONVERGENCE_REL_EPS 1.0e-30
 
typedef struct SolutionConvergenceDeterministicPass1 {
    PetscReal fluid_volume;
    PetscReal current_speed_sum;
    PetscReal current_ke_sum;
    PetscReal reference_speed_sum;
    PetscReal reference_ke_sum;
    PetscReal current_u_norm_sq;
    PetscReal delta_u_norm_sq;
    PetscReal current_pressure_sum;
    PetscReal reference_pressure_sum;
} SolutionConvergenceDeterministicPass1;
 
typedef struct SolutionConvergenceDeterministicPass2 {
    PetscReal current_pressure_norm_sq;
    PetscReal delta_pressure_norm_sq;
} SolutionConvergenceDeterministicPass2;
 
/**
 * @brief Forms a guarded relative metric for solution-convergence logging.
 *
 * This helper centralizes the divide-by-nearly-zero protection used by the
 * solution-convergence logger when turning an absolute drift into a relative
 * one. The denominator is clamped away from zero so warmup rows, quiescent
 * fields, and statistically small observables do not generate infinities.
 *
 * @param[in] numerator   Absolute quantity or drift magnitude.
 * @param[in] denominator Reference magnitude used for normalization.
 * @return Guarded relative value `numerator / max(|denominator|, eps)`.
 */
static PetscReal SolutionConvergenceSafeRelative(PetscReal numerator, PetscReal denominator)
{
    return numerator / PetscMax(PetscAbsReal(denominator), SOLUTION_CONVERGENCE_REL_EPS);
}
 
/**
 * @brief Computes instantaneous global flow observables for statistical mode.
 *
 * The statistical solution-convergence path does not compare full Eulerian
 * fields. Instead, it tracks a compact history of global observables derived
 * from the completed Eulerian state. This helper computes the current
 * volume-weighted fluid-domain mean speed and mean kinetic energy from `Ucat`.
 *
 * Only physical fluid cells contribute:
 * - solid/immersed cells are excluded using `Nvert`
 * - cells with near-zero metric Jacobian are ignored to avoid invalid volume
 *   weights
 *
 * Each MPI rank accumulates local partial sums and the routine reduces them to
 * one global pair of observables.
 *
 * @param[in]  simCtx          Simulation context owning the finest-level flow
 *                             fields.
 * @param[out] mean_speed_out  Volume-weighted domain mean of `|u|`.
 * @param[out] mean_ke_out     Volume-weighted domain mean of `0.5 |u|^2`.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode ComputeCurrentFlowObservables(SimCtx *simCtx, PetscReal *mean_speed_out, PetscReal *mean_ke_out)
{
    PetscReal local[3] = {0.0, 0.0, 0.0};
    PetscReal global[3] = {0.0, 0.0, 0.0};
    UserCtx   *user = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx || !mean_speed_out || !mean_ke_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeCurrentFlowObservables received a NULL argument.");
    }
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
        const DMDALocalInfo info = user[bi].info;
        /* Physical cells only. Index 0 and mx-1 are ghost layers on every axis - a
           boundary-condition image on a wall, a copy of the opposite cell on a
           periodic axis - and counting them biased every mean and norm here. */
        const PetscInt i_start = PetscMax(info.xs, 1), i_end = PetscMin(info.xs + info.xm, info.mx - 1);
        const PetscInt j_start = PetscMax(info.ys, 1), j_end = PetscMin(info.ys + info.ym, info.my - 1);
        const PetscInt k_start = PetscMax(info.zs, 1), k_end = PetscMin(info.zs + info.zm, info.mz - 1);
        Cmpnts           ***ucat = NULL;
        PetscReal        ***aj = NULL;
        PetscReal        ***nvert = NULL;
 
        PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
 
        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) {
                    PetscReal jac = aj[k][j][i];
                    PetscReal cell_volume = 0.0;
                    PetscReal speed = 0.0;
                    PetscReal ke = 0.0;
 
                    if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
                    if (PetscAbsReal(jac) <= 1.0e-14) continue;
 
                    cell_volume = 1.0 / jac;
                    speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
                                          ucat[k][j][i].y * ucat[k][j][i].y +
                                          ucat[k][j][i].z * ucat[k][j][i].z);
                    ke = 0.5 * speed * speed;
 
                    local[0] += cell_volume;
                    local[1] += speed * cell_volume;
                    local[2] += ke * cell_volume;
                }
            }
        }
 
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
    }
 
    PetscCallMPI(MPI_Allreduce(local, global, 3, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 
    if (global[0] <= 0.0) {
        *mean_speed_out = 0.0;
        *mean_ke_out = 0.0;
    } else {
        *mean_speed_out = global[1] / global[0];
        *mean_ke_out = global[2] / global[0];
    }
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Computes deterministic solution-drift metrics for the current step.
 *
 * This helper powers both `steady_deterministic` and
 * `periodic_deterministic` solution-convergence modes. It compares the
 * completed Eulerian state against either:
 * - the previous physical timestep (`Ucat_o`, `P_o`) for steady/transient use
 * - the stored phase-aligned snapshot ring for periodic use
 *
 * The routine performs two passes over the fluid cells:
 * 1. velocity/observable pass
 *    - computes current and reference mean speed / mean KE
 *    - computes absolute/relative velocity L2 drift
 *    - accumulates pressure means needed for gauge removal
 * 2. pressure-only pass
 *    - subtracts the volume-weighted mean pressure from both states
 *    - computes gauge-invariant pressure L2 drift
 *
 * Warmup behavior is handled here. If no valid reference exists yet, all drift
 * outputs are left at zero and `has_reference_out` is set to `PETSC_FALSE`,
 * while current observables are still reported.
 *
 * @param[in]  simCtx               Simulation context owning the current state.
 * @param[in]  periodic_mode        `PETSC_TRUE` when comparing against
 *                                  phase-aligned periodic storage.
 * @param[in]  phase_step           Active phase slot for periodic mode, or `-1`
 *                                  when unused.
 * @param[in]  samples_before       Number of solution-convergence samples
 *                                  already recorded before this timestep.
 * @param[out] has_reference_out    Whether a valid comparison state existed.
 * @param[out] u_abs_l2_out         Absolute L2 drift of Cartesian velocity.
 * @param[out] u_rel_l2_out         Relative L2 drift of Cartesian velocity.
 * @param[out] p_abs_l2_out         Gauge-invariant absolute L2 pressure drift.
 * @param[out] p_rel_l2_out         Gauge-invariant relative L2 pressure drift.
 * @param[out] mean_speed_out       Current volume-weighted mean speed.
 * @param[out] mean_speed_ref_out   Reference volume-weighted mean speed.
 * @param[out] mean_speed_abs_out   Absolute drift of mean speed.
 * @param[out] mean_speed_rel_out   Relative drift of mean speed.
 * @param[out] mean_ke_out          Current volume-weighted mean kinetic energy.
 * @param[out] mean_ke_ref_out      Reference volume-weighted mean kinetic
 *                                  energy.
 * @param[out] mean_ke_abs_out      Absolute drift of mean kinetic energy.
 * @param[out] mean_ke_rel_out      Relative drift of mean kinetic energy.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode ComputeDeterministicSolutionMetrics(SimCtx *simCtx,
                                                          PetscBool periodic_mode,
                                                          PetscInt phase_step,
                                                          PetscInt samples_before,
                                                          PetscBool *has_reference_out,
                                                          PetscReal *u_abs_l2_out,
                                                          PetscReal *u_rel_l2_out,
                                                          PetscReal *p_abs_l2_out,
                                                          PetscReal *p_rel_l2_out,
                                                          PetscReal *mean_speed_out,
                                                          PetscReal *mean_speed_ref_out,
                                                          PetscReal *mean_speed_abs_out,
                                                          PetscReal *mean_speed_rel_out,
                                                          PetscReal *mean_ke_out,
                                                          PetscReal *mean_ke_ref_out,
                                                          PetscReal *mean_ke_abs_out,
                                                          PetscReal *mean_ke_rel_out)
{
    SolutionConvergenceDeterministicPass1 local_pass1 = {0};
    SolutionConvergenceDeterministicPass1 global_pass1 = {0};
    SolutionConvergenceDeterministicPass2 local_pass2 = {0};
    SolutionConvergenceDeterministicPass2 global_pass2 = {0};
    PetscReal current_pressure_mean = 0.0;
    PetscReal reference_pressure_mean = 0.0;
    UserCtx   *user = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx || !has_reference_out || !u_abs_l2_out || !u_rel_l2_out || !p_abs_l2_out || !p_rel_l2_out ||
        !mean_speed_out || !mean_speed_ref_out || !mean_speed_abs_out || !mean_speed_rel_out ||
        !mean_ke_out || !mean_ke_ref_out || !mean_ke_abs_out || !mean_ke_rel_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeDeterministicSolutionMetrics received a NULL output pointer.");
    }
 
    *has_reference_out = periodic_mode
        ? (PetscBool)(simCtx->solutionConvergencePeriodSteps > 0 &&
                      phase_step >= 0 &&
                      phase_step < simCtx->solutionConvergencePeriodSteps &&
                      samples_before >= simCtx->solutionConvergencePeriodSteps)
        : (PetscBool)(samples_before > 0);
 
    *u_abs_l2_out = 0.0;
    *u_rel_l2_out = 0.0;
    *p_abs_l2_out = 0.0;
    *p_rel_l2_out = 0.0;
    *mean_speed_out = 0.0;
    *mean_speed_ref_out = 0.0;
    *mean_speed_abs_out = 0.0;
    *mean_speed_rel_out = 0.0;
    *mean_ke_out = 0.0;
    *mean_ke_ref_out = 0.0;
    *mean_ke_abs_out = 0.0;
    *mean_ke_rel_out = 0.0;
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
        const DMDALocalInfo info = user[bi].info;
        /* Physical cells only. Index 0 and mx-1 are ghost layers on every axis - a
           boundary-condition image on a wall, a copy of the opposite cell on a
           periodic axis - and counting them biased every mean and norm here. */
        const PetscInt i_start = PetscMax(info.xs, 1), i_end = PetscMin(info.xs + info.xm, info.mx - 1);
        const PetscInt j_start = PetscMax(info.ys, 1), j_end = PetscMin(info.ys + info.ym, info.my - 1);
        const PetscInt k_start = PetscMax(info.zs, 1), k_end = PetscMin(info.zs + info.zm, info.mz - 1);
        Cmpnts           ***ucat = NULL;
        Cmpnts           ***ucat_ref = NULL;
        PetscReal        ***pressure = NULL;
        PetscReal        ***pressure_ref = NULL;
        PetscReal        ***aj = NULL;
        PetscReal        ***nvert = NULL;
        Vec                 ucat_reference_vec = NULL;
        Vec                 pressure_reference_vec = NULL;
 
        if (*has_reference_out) {
            if (periodic_mode) {
                ucat_reference_vec = user[bi].solutionConvergencePeriodicUcatRef[phase_step];
                pressure_reference_vec = user[bi].solutionConvergencePeriodicPRef[phase_step];
            } else {
                ucat_reference_vec = user[bi].Ucat_o;
                pressure_reference_vec = user[bi].P_o;
            }
        }
 
        PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
        if (*has_reference_out) {
            PetscCall(DMDAVecGetArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
            PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
        }
 
        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) {
                    PetscReal jac = aj[k][j][i];
                    PetscReal cell_volume = 0.0;
                    PetscReal speed = 0.0;
                    PetscReal ke = 0.0;
 
                    if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
                    if (PetscAbsReal(jac) <= 1.0e-14) continue;
 
                    cell_volume = 1.0 / jac;
                    speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
                                          ucat[k][j][i].y * ucat[k][j][i].y +
                                          ucat[k][j][i].z * ucat[k][j][i].z);
                    ke = 0.5 * speed * speed;
 
                    local_pass1.fluid_volume += cell_volume;
                    local_pass1.current_speed_sum += speed * cell_volume;
                    local_pass1.current_ke_sum += ke * cell_volume;
                    local_pass1.current_u_norm_sq += (ucat[k][j][i].x * ucat[k][j][i].x +
                                                      ucat[k][j][i].y * ucat[k][j][i].y +
                                                      ucat[k][j][i].z * ucat[k][j][i].z) * cell_volume;
 
                    if (*has_reference_out) {
                        PetscReal ref_speed = PetscSqrtReal(ucat_ref[k][j][i].x * ucat_ref[k][j][i].x +
                                                            ucat_ref[k][j][i].y * ucat_ref[k][j][i].y +
                                                            ucat_ref[k][j][i].z * ucat_ref[k][j][i].z);
                        PetscReal ref_ke = 0.5 * ref_speed * ref_speed;
                        PetscReal dux = ucat[k][j][i].x - ucat_ref[k][j][i].x;
                        PetscReal duy = ucat[k][j][i].y - ucat_ref[k][j][i].y;
                        PetscReal duz = ucat[k][j][i].z - ucat_ref[k][j][i].z;
 
                        local_pass1.reference_speed_sum += ref_speed * cell_volume;
                        local_pass1.reference_ke_sum += ref_ke * cell_volume;
                        local_pass1.delta_u_norm_sq += (dux * dux + duy * duy + duz * duz) * cell_volume;
                        local_pass1.current_pressure_sum += pressure[k][j][i] * cell_volume;
                        local_pass1.reference_pressure_sum += pressure_ref[k][j][i] * cell_volume;
                    }
                }
            }
        }
 
        if (*has_reference_out) {
            PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
            PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
        }
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
    }
 
    PetscCallMPI(MPI_Allreduce(&local_pass1, &global_pass1,
                               sizeof(SolutionConvergenceDeterministicPass1) / sizeof(PetscReal),
                               MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 
    if (global_pass1.fluid_volume <= 0.0) PetscFunctionReturn(0);
 
    *mean_speed_out = global_pass1.current_speed_sum / global_pass1.fluid_volume;
    *mean_ke_out = global_pass1.current_ke_sum / global_pass1.fluid_volume;
 
    if (!*has_reference_out) PetscFunctionReturn(0);
 
    *mean_speed_ref_out = global_pass1.reference_speed_sum / global_pass1.fluid_volume;
    *mean_speed_abs_out = PetscAbsReal(*mean_speed_out - *mean_speed_ref_out);
    *mean_speed_rel_out = SolutionConvergenceSafeRelative(*mean_speed_abs_out, *mean_speed_out);
    *mean_ke_ref_out = global_pass1.reference_ke_sum / global_pass1.fluid_volume;
    *mean_ke_abs_out = PetscAbsReal(*mean_ke_out - *mean_ke_ref_out);
    *mean_ke_rel_out = SolutionConvergenceSafeRelative(*mean_ke_abs_out, *mean_ke_out);
 
    current_pressure_mean = global_pass1.current_pressure_sum / global_pass1.fluid_volume;
    reference_pressure_mean = global_pass1.reference_pressure_sum / global_pass1.fluid_volume;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
        const DMDALocalInfo info = user[bi].info;
        /* Physical cells only. Index 0 and mx-1 are ghost layers on every axis - a
           boundary-condition image on a wall, a copy of the opposite cell on a
           periodic axis - and counting them biased every mean and norm here. */
        const PetscInt i_start = PetscMax(info.xs, 1), i_end = PetscMin(info.xs + info.xm, info.mx - 1);
        const PetscInt j_start = PetscMax(info.ys, 1), j_end = PetscMin(info.ys + info.ym, info.my - 1);
        const PetscInt k_start = PetscMax(info.zs, 1), k_end = PetscMin(info.zs + info.zm, info.mz - 1);
        PetscReal        ***pressure = NULL;
        PetscReal        ***pressure_ref = NULL;
        PetscReal        ***aj = NULL;
        PetscReal        ***nvert = NULL;
        Vec                 pressure_reference_vec = periodic_mode ? user[bi].solutionConvergencePeriodicPRef[phase_step] : user[bi].P_o;
 
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
 
        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) {
                    PetscReal jac = aj[k][j][i];
                    PetscReal cell_volume = 0.0;
                    PetscReal current_pressure = 0.0;
                    PetscReal reference_pressure = 0.0;
                    PetscReal delta_pressure = 0.0;
 
                    if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
                    if (PetscAbsReal(jac) <= 1.0e-14) continue;
 
                    cell_volume = 1.0 / jac;
                    current_pressure = pressure[k][j][i] - current_pressure_mean;
                    reference_pressure = pressure_ref[k][j][i] - reference_pressure_mean;
                    delta_pressure = current_pressure - reference_pressure;
 
                    local_pass2.current_pressure_norm_sq += current_pressure * current_pressure * cell_volume;
                    local_pass2.delta_pressure_norm_sq += delta_pressure * delta_pressure * cell_volume;
                }
            }
        }
 
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
        PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
    }
 
    PetscCallMPI(MPI_Allreduce(&local_pass2, &global_pass2,
                               sizeof(SolutionConvergenceDeterministicPass2) / sizeof(PetscReal),
                               MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 
    *u_abs_l2_out = PetscSqrtReal(global_pass1.delta_u_norm_sq);
    *u_rel_l2_out = SolutionConvergenceSafeRelative(*u_abs_l2_out, PetscSqrtReal(global_pass1.current_u_norm_sq));
    *p_abs_l2_out = PetscSqrtReal(global_pass2.delta_pressure_norm_sq);
    *p_rel_l2_out = SolutionConvergenceSafeRelative(*p_abs_l2_out, PetscSqrtReal(global_pass2.current_pressure_norm_sq));
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Reads one sample from the statistical ring buffer by age.
 *
 * The statistical logger stores scalar observables in a compact circular
 * buffer. This helper interprets the buffer using `samples_available` as the
 * logical end of the history and returns the entry `offset_from_latest` steps
 * back from the newest stored sample.
 *
 * Out-of-range requests return zero so warmup handling can remain simple and
 * deterministic.
 *
 * @param[in] history              Ring-buffer storage array.
 * @param[in] capacity             Total ring-buffer capacity.
 * @param[in] samples_available    Number of logical samples available to read.
 * @param[in] offset_from_latest   `0` means newest sample, `1` previous sample,
 *                                 and so on.
 * @return Requested historical sample, or `0.0` if unavailable.
 */
static PetscReal SolutionConvergenceHistoryGet(const PetscReal *history,
                                               PetscInt capacity,
                                               PetscInt samples_available,
                                               PetscInt offset_from_latest)
{
    PetscInt count = 0;
    PetscInt index = 0;
 
    if (!history || capacity <= 0 || samples_available <= 0 || offset_from_latest < 0) {
        return 0.0;
    }
 
    count = PetscMin(samples_available, capacity);
    if (offset_from_latest >= count) return 0.0;
 
    index = (samples_available - 1 - offset_from_latest) % capacity;
    if (index < 0) index += capacity;
    return history[index];
}
 
/**
 * @brief Appends one timestep's scalar observables to the statistical history.
 *
 * Statistical solution-convergence compares adjacent windows of scalar
 * observables rather than full fields. This helper writes the current
 * `mean_speed` and `mean_ke` into the rolling history arrays using the current
 * sample count to choose the circular-buffer slot.
 *
 * @param[in,out] simCtx         Simulation context owning the history arrays.
 * @param[in]     samples_before Number of samples present before appending the
 *                               current timestep.
 * @param[in]     mean_speed     Current timestep mean-speed observable.
 * @param[in]     mean_ke        Current timestep mean-KE observable.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode AppendStatisticalObservableSample(SimCtx *simCtx,
                                                        PetscInt samples_before,
                                                        PetscReal mean_speed,
                                                        PetscReal mean_ke)
{
    PetscInt history_capacity = 0;
    PetscInt slot = 0;
 
    PetscFunctionBeginUser;
    if (!simCtx || !simCtx->solutionConvergenceMeanSpeedHistory || !simCtx->solutionConvergenceMeanKEHistory) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Statistical solution-convergence history is not allocated.");
    }
 
    history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
    if (history_capacity <= 0) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Statistical solution-convergence history capacity must be positive.");
    }
 
    slot = samples_before % history_capacity;
    simCtx->solutionConvergenceMeanSpeedHistory[slot] = mean_speed;
    simCtx->solutionConvergenceMeanKEHistory[slot] = mean_ke;
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Computes adjacent-window drift metrics for statistical steady mode.
 *
 * Once enough samples have been accumulated, this helper forms:
 * - the current window over the most recent `window_steps` samples
 * - the previous adjacent window over the preceding `window_steps` samples
 *
 * From those windows it computes means, RMS values, and absolute/relative
 * drift for both tracked observables (`mean_speed` and `mean_ke`). When the
 * history is still warming up:
 * - fewer than `window_steps` samples: no window metrics are available
 * - between `window_steps` and `2*window_steps - 1` samples: current-window
 *   metrics are available, but no reference window exists yet
 *
 * @param[in]  simCtx                           Simulation context owning the
 *                                             statistical history.
 * @param[in]  samples_available               Number of samples available after
 *                                             appending the current timestep.
 * @param[out] has_reference_out               Whether both adjacent windows
 *                                             exist and drift metrics are
 *                                             meaningful.
 * @param[out] mean_speed_window_out           Mean speed over the current
 *                                             window.
 * @param[out] mean_speed_window_prev_out      Mean speed over the previous
 *                                             window.
 * @param[out] mean_speed_window_abs_out       Absolute drift between current
 *                                             and previous window means.
 * @param[out] mean_speed_window_rel_out       Relative drift between current
 *                                             and previous window means.
 * @param[out] mean_speed_rms_window_out       RMS of mean-speed samples in the
 *                                             current window.
 * @param[out] mean_speed_rms_window_prev_out  RMS of mean-speed samples in the
 *                                             previous window.
 * @param[out] mean_speed_rms_window_abs_out   Absolute drift between window RMS
 *                                             values.
 * @param[out] mean_speed_rms_window_rel_out   Relative drift between window RMS
 *                                             values.
 * @param[out] mean_ke_window_out              Mean kinetic energy over the
 *                                             current window.
 * @param[out] mean_ke_window_prev_out         Mean kinetic energy over the
 *                                             previous window.
 * @param[out] mean_ke_window_abs_out          Absolute drift between current
 *                                             and previous KE-window means.
 * @param[out] mean_ke_window_rel_out          Relative drift between current
 *                                             and previous KE-window means.
 * @param[out] mean_ke_rms_window_out          RMS of mean-KE samples in the
 *                                             current window.
 * @param[out] mean_ke_rms_window_prev_out     RMS of mean-KE samples in the
 *                                             previous window.
 * @param[out] mean_ke_rms_window_abs_out      Absolute drift between KE-window
 *                                             RMS values.
 * @param[out] mean_ke_rms_window_rel_out      Relative drift between KE-window
 *                                             RMS values.
 * @return PetscErrorCode 0 on success.
 */
static PetscErrorCode ComputeStatisticalWindowMetrics(const SimCtx *simCtx,
                                                      PetscInt samples_available,
                                                      PetscBool *has_reference_out,
                                                      PetscReal *mean_speed_window_out,
                                                      PetscReal *mean_speed_window_prev_out,
                                                      PetscReal *mean_speed_window_abs_out,
                                                      PetscReal *mean_speed_window_rel_out,
                                                      PetscReal *mean_speed_rms_window_out,
                                                      PetscReal *mean_speed_rms_window_prev_out,
                                                      PetscReal *mean_speed_rms_window_abs_out,
                                                      PetscReal *mean_speed_rms_window_rel_out,
                                                      PetscReal *mean_ke_window_out,
                                                      PetscReal *mean_ke_window_prev_out,
                                                      PetscReal *mean_ke_window_abs_out,
                                                      PetscReal *mean_ke_window_rel_out,
                                                      PetscReal *mean_ke_rms_window_out,
                                                      PetscReal *mean_ke_rms_window_prev_out,
                                                      PetscReal *mean_ke_rms_window_abs_out,
                                                      PetscReal *mean_ke_rms_window_rel_out)
{
    PetscInt  w = 0;
    PetscInt  history_capacity = 0;
    PetscReal speed_sum = 0.0;
    PetscReal speed_sum_sq = 0.0;
    PetscReal speed_prev_sum = 0.0;
    PetscReal speed_prev_sum_sq = 0.0;
    PetscReal ke_sum = 0.0;
    PetscReal ke_sum_sq = 0.0;
    PetscReal ke_prev_sum = 0.0;
    PetscReal ke_prev_sum_sq = 0.0;
 
    PetscFunctionBeginUser;
    if (!simCtx || !has_reference_out || !mean_speed_window_out || !mean_speed_window_prev_out ||
        !mean_speed_window_abs_out || !mean_speed_window_rel_out || !mean_speed_rms_window_out ||
        !mean_speed_rms_window_prev_out || !mean_speed_rms_window_abs_out || !mean_speed_rms_window_rel_out ||
        !mean_ke_window_out || !mean_ke_window_prev_out || !mean_ke_window_abs_out || !mean_ke_window_rel_out ||
        !mean_ke_rms_window_out || !mean_ke_rms_window_prev_out || !mean_ke_rms_window_abs_out ||
        !mean_ke_rms_window_rel_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeStatisticalWindowMetrics received a NULL output pointer.");
    }
 
    *has_reference_out = PETSC_FALSE;
    *mean_speed_window_out = 0.0;
    *mean_speed_window_prev_out = 0.0;
    *mean_speed_window_abs_out = 0.0;
    *mean_speed_window_rel_out = 0.0;
    *mean_speed_rms_window_out = 0.0;
    *mean_speed_rms_window_prev_out = 0.0;
    *mean_speed_rms_window_abs_out = 0.0;
    *mean_speed_rms_window_rel_out = 0.0;
    *mean_ke_window_out = 0.0;
    *mean_ke_window_prev_out = 0.0;
    *mean_ke_window_abs_out = 0.0;
    *mean_ke_window_rel_out = 0.0;
    *mean_ke_rms_window_out = 0.0;
    *mean_ke_rms_window_prev_out = 0.0;
    *mean_ke_rms_window_abs_out = 0.0;
    *mean_ke_rms_window_rel_out = 0.0;
 
    w = simCtx->solutionConvergenceWindowSteps;
    history_capacity = 2 * w;
    if (w <= 0 || samples_available < w) PetscFunctionReturn(0);
 
    for (PetscInt idx = 0; idx < w; ++idx) {
        PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
                                                              history_capacity,
                                                              samples_available,
                                                              idx);
        PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
                                                           history_capacity,
                                                           samples_available,
                                                           idx);
        speed_sum += speed_value;
        speed_sum_sq += speed_value * speed_value;
        ke_sum += ke_value;
        ke_sum_sq += ke_value * ke_value;
    }
 
    *mean_speed_window_out = speed_sum / (PetscReal)w;
    *mean_speed_rms_window_out = PetscSqrtReal(PetscMax(0.0, speed_sum_sq / (PetscReal)w -
                                                             (*mean_speed_window_out) * (*mean_speed_window_out)));
    *mean_ke_window_out = ke_sum / (PetscReal)w;
    *mean_ke_rms_window_out = PetscSqrtReal(PetscMax(0.0, ke_sum_sq / (PetscReal)w -
                                                          (*mean_ke_window_out) * (*mean_ke_window_out)));
 
    if (samples_available < 2 * w) PetscFunctionReturn(0);
 
    for (PetscInt idx = w; idx < 2 * w; ++idx) {
        PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
                                                              history_capacity,
                                                              samples_available,
                                                              idx);
        PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
                                                           history_capacity,
                                                           samples_available,
                                                           idx);
        speed_prev_sum += speed_value;
        speed_prev_sum_sq += speed_value * speed_value;
        ke_prev_sum += ke_value;
        ke_prev_sum_sq += ke_value * ke_value;
    }
 
    *has_reference_out = PETSC_TRUE;
    *mean_speed_window_prev_out = speed_prev_sum / (PetscReal)w;
    *mean_speed_window_abs_out = PetscAbsReal(*mean_speed_window_out - *mean_speed_window_prev_out);
    *mean_speed_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_window_abs_out, *mean_speed_window_out);
    *mean_speed_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, speed_prev_sum_sq / (PetscReal)w -
                                                                   (*mean_speed_window_prev_out) * (*mean_speed_window_prev_out)));
    *mean_speed_rms_window_abs_out = PetscAbsReal(*mean_speed_rms_window_out - *mean_speed_rms_window_prev_out);
    *mean_speed_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_rms_window_abs_out, *mean_speed_rms_window_out);
 
    *mean_ke_window_prev_out = ke_prev_sum / (PetscReal)w;
    *mean_ke_window_abs_out = PetscAbsReal(*mean_ke_window_out - *mean_ke_window_prev_out);
    *mean_ke_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_window_abs_out, *mean_ke_window_out);
    *mean_ke_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, ke_prev_sum_sq / (PetscReal)w -
                                                                (*mean_ke_window_prev_out) * (*mean_ke_window_prev_out)));
    *mean_ke_rms_window_abs_out = PetscAbsReal(*mean_ke_rms_window_out - *mean_ke_rms_window_prev_out);
    *mean_ke_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_rms_window_abs_out, *mean_ke_rms_window_out);
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Maps the internal solution-convergence mode enum to its log label.
 *
 * The logger writes a human-readable mode string into the
 * `solution_convergence.log` banner and `mode` column. This helper keeps the
 * formatting centralized so the file output stays consistent with the accepted
 * configuration names.
 *
 * @param[in] mode Internal solution-convergence mode selector.
 * @return Lowercase string label written to the log output.
 */
static const char *SolutionConvergenceModeToString(SolutionConvergenceMode mode)
{
    switch (mode) {
        case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC: return "steady_deterministic";
        case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC: return "periodic_deterministic";
        case SOLUTION_CONVERGENCE_STATISTICAL_STEADY: return "statistical_steady";
        case SOLUTION_CONVERGENCE_TRANSIENT: return "transient";
        default: return "unknown";
    }
}
 
/**
 * @brief Implementation of \ref LOG_SOLUTION_CONVERGENCE().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_SOLUTION_CONVERGENCE()
 */
PetscErrorCode LOG_SOLUTION_CONVERGENCE(SimCtx *simCtx)
{
    PetscMPIInt rank = 0;
    PetscBool   has_reference = PETSC_FALSE;
    PetscInt    phase_step = -1;
    PetscInt    samples_before = 0;
    PetscReal   u_abs_l2 = 0.0, u_rel_l2 = 0.0, p_abs_l2 = 0.0, p_rel_l2 = 0.0;
    PetscReal   mean_speed = 0.0, mean_speed_reference = 0.0, mean_speed_abs_drift = 0.0, mean_speed_rel_drift = 0.0;
    PetscReal   mean_ke = 0.0, mean_ke_reference = 0.0, mean_ke_abs_drift = 0.0, mean_ke_rel_drift = 0.0;
    PetscReal   mean_speed_window = 0.0, mean_speed_window_prev = 0.0, mean_speed_window_abs_drift = 0.0, mean_speed_window_rel_drift = 0.0;
    PetscReal   mean_speed_rms_window = 0.0, mean_speed_rms_window_prev = 0.0, mean_speed_rms_window_abs_drift = 0.0, mean_speed_rms_window_rel_drift = 0.0;
    PetscReal   mean_ke_window = 0.0, mean_ke_window_prev = 0.0, mean_ke_window_abs_drift = 0.0, mean_ke_window_rel_drift = 0.0;
    PetscReal   mean_ke_rms_window = 0.0, mean_ke_rms_window_prev = 0.0, mean_ke_rms_window_abs_drift = 0.0, mean_ke_rms_window_rel_drift = 0.0;
 
    PetscFunctionBeginUser;
    if (!simCtx) PetscFunctionReturn(0);
    if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
    if (!simCtx->solutionConvergenceEnabled) PetscFunctionReturn(0);
 
    samples_before = simCtx->solutionConvergenceSamplesRecorded;
 
    switch (simCtx->solutionConvergenceMode) {
        case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
        case SOLUTION_CONVERGENCE_TRANSIENT:
            PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_FALSE, -1, samples_before,
                                                          &has_reference,
                                                          &u_abs_l2, &u_rel_l2,
                                                          &p_abs_l2, &p_rel_l2,
                                                          &mean_speed, &mean_speed_reference,
                                                          &mean_speed_abs_drift, &mean_speed_rel_drift,
                                                          &mean_ke, &mean_ke_reference,
                                                          &mean_ke_abs_drift, &mean_ke_rel_drift));
            break;
        case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
            phase_step = simCtx->solutionConvergencePeriodSteps > 0 ? (simCtx->step % simCtx->solutionConvergencePeriodSteps) : -1;
            PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_TRUE, phase_step, samples_before,
                                                          &has_reference,
                                                          &u_abs_l2, &u_rel_l2,
                                                          &p_abs_l2, &p_rel_l2,
                                                          &mean_speed, &mean_speed_reference,
                                                          &mean_speed_abs_drift, &mean_speed_rel_drift,
                                                          &mean_ke, &mean_ke_reference,
                                                          &mean_ke_abs_drift, &mean_ke_rel_drift));
            break;
        case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
            PetscCall(ComputeCurrentFlowObservables(simCtx, &mean_speed, &mean_ke));
            PetscCall(AppendStatisticalObservableSample(simCtx, samples_before, mean_speed, mean_ke));
            PetscCall(ComputeStatisticalWindowMetrics(simCtx, samples_before + 1,
                                                      &has_reference,
                                                      &mean_speed_window, &mean_speed_window_prev,
                                                      &mean_speed_window_abs_drift, &mean_speed_window_rel_drift,
                                                      &mean_speed_rms_window, &mean_speed_rms_window_prev,
                                                      &mean_speed_rms_window_abs_drift, &mean_speed_rms_window_rel_drift,
                                                      &mean_ke_window, &mean_ke_window_prev,
                                                      &mean_ke_window_abs_drift, &mean_ke_window_rel_drift,
                                                      &mean_ke_rms_window, &mean_ke_rms_window_prev,
                                                      &mean_ke_rms_window_abs_drift, &mean_ke_rms_window_rel_drift));
            break;
        default:
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown solution convergence mode %d.", (int)simCtx->solutionConvergenceMode);
    }
 
    PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    if (rank == 0) {
        char log_path[PETSC_MAX_PATH_LEN + 32];
        FILE *f = NULL;
        const char *mode_str = SolutionConvergenceModeToString(simCtx->solutionConvergenceMode);
 
        PetscCall(PetscSNPrintf(log_path, sizeof(log_path), "%s/solution_convergence.log", simCtx->log_dir));
        f = fopen(log_path, "a");
        if (!f) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open solution convergence log: %s", log_path);
        }
 
        if (ftell(f) == 0) {
            switch (simCtx->solutionConvergenceMode) {
                case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
                case SOLUTION_CONVERGENCE_TRANSIENT:
                    fprintf(f, "==================== Solution Convergence Log [mode: %s] ====================\n", mode_str);
                    /* 16 columns; header width = 314 chars */
                    fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
                            "step", "time", "mode", "ref",
                            "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
                            "mean_speed", "spd_ref", "spd_abs", "spd_rel",
                            "mean_ke", "ke_ref", "ke_abs", "ke_rel");
                    fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
                    break;
                case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
                    fprintf(f, "==================== Solution Convergence Log [mode: %s | period_steps: %d] ====================\n",
                            mode_str, (int)simCtx->solutionConvergencePeriodSteps);
                    /* 18 columns; header width = 330 chars */
                    fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
                            "step", "time", "mode", "ref", "ph", "per",
                            "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
                            "mean_speed", "spd_ref", "spd_abs", "spd_rel",
                            "mean_ke", "ke_ref", "ke_abs", "ke_rel");
                    fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
                    break;
                case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
                    fprintf(f, "==================== Solution Convergence Log [mode: %s | window_steps: %d] ====================\n",
                            mode_str, (int)simCtx->solutionConvergenceWindowSteps);
                    /* 21 columns; header width = 406 chars */
                    fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
                            "step", "time", "mode", "ref", "win",
                            "mean_speed", "mean_ke",
                            "spd_win", "spd_win_prev", "spd_win_abs", "spd_win_rel",
                            "spd_rms_win", "spd_rms_abs", "spd_rms_rel",
                            "ke_win", "ke_win_prev", "ke_win_abs", "ke_win_rel",
                            "ke_rms_win", "ke_rms_abs", "ke_rms_rel");
                    fprintf(f, "------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
                    break;
                default: break;
            }
        }
        if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
            fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
        }
 
        switch (simCtx->solutionConvergenceMode) {
            case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
            case SOLUTION_CONVERGENCE_TRANSIENT:
                fprintf(f,
                        "%-10d | %-18.10e | %-22s | %-3d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
                        (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
                        (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
                        (double)mean_speed, (double)mean_speed_reference,
                        (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
                        (double)mean_ke, (double)mean_ke_reference,
                        (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
                break;
            case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
                fprintf(f,
                        "%-10d | %-18.10e | %-22s | %-3d | %-5d | %-5d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
                        (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
                        (int)phase_step, (int)simCtx->solutionConvergencePeriodSteps,
                        (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
                        (double)mean_speed, (double)mean_speed_reference,
                        (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
                        (double)mean_ke, (double)mean_ke_reference,
                        (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
                break;
            case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
                fprintf(f,
                        "%-10d | %-18.10e | %-22s | %-3d | %-5d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
                        (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
                        (int)simCtx->solutionConvergenceWindowSteps,
                        (double)mean_speed, (double)mean_ke,
                        (double)mean_speed_window, (double)mean_speed_window_prev,
                        (double)mean_speed_window_abs_drift, (double)mean_speed_window_rel_drift,
                        (double)mean_speed_rms_window,
                        (double)mean_speed_rms_window_abs_drift, (double)mean_speed_rms_window_rel_drift,
                        (double)mean_ke_window, (double)mean_ke_window_prev,
                        (double)mean_ke_window_abs_drift, (double)mean_ke_window_rel_drift,
                        (double)mean_ke_rms_window,
                        (double)mean_ke_rms_window_abs_drift, (double)mean_ke_rms_window_rel_drift);
                break;
            default: break;
        }
        fclose(f);
    }
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
        phase_step >= 0 && phase_step < simCtx->solutionConvergencePeriodSteps) {
        UserCtx *user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            PetscCall(VecCopy(user[bi].Ucat, user[bi].solutionConvergencePeriodicUcatRef[phase_step]));
            PetscCall(VecCopy(user[bi].P, user[bi].solutionConvergencePeriodicPRef[phase_step]));
        }
    }
 
    simCtx->solutionConvergenceSamplesRecorded = samples_before + 1;
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Logs continuity metrics for a single block to a file.
 *
 * This function should be called for each block, once per timestep. It opens a
 * central log file in append mode. To ensure the header is written only once,
 * it checks if it is processing block 0 on the simulation's start step.
 *
 * @param user A pointer to the UserCtx for the specific block whose metrics
 *             are to be logged. The function accesses both global (SimCtx)
 *             and local (user->...) data.
 * @return     PetscErrorCode 0 on success.
 */
#undef __FUNCT__
#define __FUNCT__ "LOG_CONTINUITY_METRICS"
/**
 * @brief Implementation of \ref LOG_CONTINUITY_METRICS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_CONTINUITY_METRICS()
 */
 
PetscErrorCode LOG_CONTINUITY_METRICS(UserCtx *user)
{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    SimCtx         *simCtx = user->simCtx;      // Get the shared SimCtx
    const PetscInt bi = user->_this;          // Get this block's specific ID
    const PetscInt ti = simCtx->step;         // Get the current timestep
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // Only rank 0 performs file I/O.
    if (!rank) {
        FILE *f;
        char filen[PETSC_MAX_PATH_LEN + 64];
        ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Continuity_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
 
        // Open the log file in append mode.
        f = fopen(filen, "a");
        if (!f) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
        }
 
        // Write a header only when the file is empty and it's the first block (bi=0).
        // Using ftell() instead of step comparison ensures correctness across continuations.
        if (ftell(f) == 0 && bi == 0) {
            PetscFPrintf(PETSC_COMM_SELF, f, "%-10s | %-6s | %-18s | %-30s | %-24s | %-18s | %-18s | %-18s\n",
                         "Timestep", "Block", "Max Divergence", "Max Divergence Location ([k][j][i]=idx)", "Poisson Source Imbalance","Total Flux In", "Total Flux Out", "Net Flux");
            PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------------------------------\n");
        }
        if (simCtx->continueMode && ti == simCtx->StartStep + 1 && bi == 0) {
            PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
        }
 
        // Prepare the data strings and values for the current block.
        PetscReal net_flux = simCtx->FluxInSum - simCtx->FluxOutSum;
        char location_str[64];
        sprintf(location_str, "([%d][%d][%d] = %d)", (int)simCtx->MaxDivz, (int)simCtx->MaxDivy, (int)simCtx->MaxDivx, (int)simCtx->MaxDivFlatArg);
 
        // Write the formatted line for the current block.
        PetscFPrintf(PETSC_COMM_SELF, f, "%-10d | %-6d | %-18.10e | %-39s | %-24.10e | %-18.10e | %-18.10e | %-18.10e\n",
                     (int)ti,
                     (int)bi,
                     (double)simCtx->MaxDiv,
                     location_str,
                     (double)simCtx->poissonSourceImbalance,
             (double)simCtx->FluxInSum,
                     (double)simCtx->FluxOutSum,
                     (double)net_flux);
 
        fclose(f);
    }
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref ParticleLocationStatusToString().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ParticleLocationStatusToString()
 */
const char* ParticleLocationStatusToString(ParticleLocationStatus level) 
{
    switch (level) {
        case NEEDS_LOCATION:        return "NEEDS_LOCATION";
        case ACTIVE_AND_LOCATED: return "ACTIVE_AND_LOCATED";
        case MIGRATING_OUT:     return "MIGRATING_OUT";
        case LOST:               return "LOST";
        case UNINITIALIZED:      return "UNINITIALIZED";
        default:            return "UNKNOWN_LEVEL";
    }
}
 
///////// Profiling System /////////
 
// Data structure to hold profiling info for one function
typedef struct {
    const char *name;
    double     total_time;
    double     current_step_time;
    long long  total_call_count;
    long long  current_step_call_count;
    double     start_time; // Timer for the current call
    PetscBool  always_log;
} ProfiledFunction;
 
// Global registry for all profiled functions
static ProfiledFunction *g_profiler_registry = NULL;
static PetscInt g_profiler_count = 0;
static PetscInt g_profiler_capacity = 0;
 
// Internal helper to find a function in the registry or create it
/**
 * @brief Find a profiling record by name or allocate and register a new record.
 */
static PetscErrorCode _FindOrCreateEntry(const char *func_name, PetscInt *idx)
{
    PetscFunctionBeginUser;
    // Search for existing entry
    for (PetscInt i = 0; i < g_profiler_count; ++i) {
        if (strcmp(g_profiler_registry[i].name, func_name) == 0) {
            *idx = i;
            PetscFunctionReturn(0);
        }
    }
 
    // Not found, create a new entry
    if (g_profiler_count >= g_profiler_capacity) {
        PetscInt new_capacity = g_profiler_capacity == 0 ? 16 : g_profiler_capacity * 2;
        PetscErrorCode ierr = PetscRealloc(sizeof(ProfiledFunction) * new_capacity, &g_profiler_registry); CHKERRQ(ierr);
        g_profiler_capacity = new_capacity;
    }
 
    *idx = g_profiler_count;
    g_profiler_registry[*idx].name = func_name;
    g_profiler_registry[*idx].total_time = 0.0;
    g_profiler_registry[*idx].current_step_time = 0.0;
    g_profiler_registry[*idx].total_call_count = 0;
    g_profiler_registry[*idx].current_step_call_count = 0;
    g_profiler_registry[*idx].start_time = 0.0;
    g_profiler_registry[*idx].always_log = PETSC_FALSE;
    g_profiler_count++;
 
    PetscFunctionReturn(0);
}
 
// --- Public API Implementation ---
/**
 * @brief Internal helper implementation: `ProfilingInitialize()`.
 * @details Local to this translation unit.
 */
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
{
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingInitialize");
 
    // Iterate through the list of critical functions provided in SimCtx
    for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; ++i) {
        PetscInt idx;
        const char *func_name = simCtx->profilingSelectedFuncs[i];
        PetscErrorCode ierr = _FindOrCreateEntry(func_name, &idx); CHKERRQ(ierr);
        g_profiler_registry[idx].always_log = PETSC_TRUE;
 
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Marked '%s' as a critical function for profiling.\n", func_name);
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref _ProfilingStart().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see _ProfilingStart()
 */
 
void _ProfilingStart(const char *func_name)
{
    PetscInt idx;
    if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
    PetscTime(&g_profiler_registry[idx].start_time);
}
 
/**
 * @brief Implementation of \ref _ProfilingEnd().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see _ProfilingEnd()
 */
 
void _ProfilingEnd(const char *func_name)
{
    double end_time;
    PetscTime(&end_time);
 
    PetscInt idx;
    if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
 
    double elapsed = end_time - g_profiler_registry[idx].start_time;
    g_profiler_registry[idx].total_time += elapsed;
    g_profiler_registry[idx].current_step_time += elapsed;
    g_profiler_registry[idx].total_call_count++;
    g_profiler_registry[idx].current_step_call_count++;
}
 
/**
 * @brief Implementation of \ref ProfilingResetTimestepCounters().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ProfilingResetTimestepCounters()
 */
 
PetscErrorCode ProfilingResetTimestepCounters(void)
{
    PetscFunctionBeginUser;
    for (PetscInt i = 0; i < g_profiler_count; ++i) {
        g_profiler_registry[i].current_step_time = 0.0;
        g_profiler_registry[i].current_step_call_count = 0;
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref ProfilingLogTimestepSummary().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ProfilingLogTimestepSummary()
 */
 
PetscErrorCode ProfilingLogTimestepSummary(SimCtx *simCtx, PetscInt step)
{
    PetscBool should_write = PETSC_FALSE;
    FILE *f = NULL;
    char filen[(2 * PETSC_MAX_PATH_LEN) + 16];
 
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingLogTimestepSummary");
 
    if (strcmp(simCtx->profilingTimestepMode, "off") == 0) {
        for (PetscInt i = 0; i < g_profiler_count; ++i) {
            g_profiler_registry[i].current_step_time = 0.0;
            g_profiler_registry[i].current_step_call_count = 0;
        }
        PetscFunctionReturn(0);
    }
 
    for (PetscInt i = 0; i < g_profiler_count; ++i) {
        if (g_profiler_registry[i].current_step_call_count <= 0) {
            continue;
        }
        if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
            should_write = PETSC_TRUE;
            break;
        }
    }
 
    /* A selected name that matches no instrumented function records nothing, and a
       list made only of such names writes no file at all. Say so once, on the first
       step, rather than leave an absent file to be discovered after the run. */
    if (step == simCtx->StartStep + 1 && strcmp(simCtx->profilingTimestepMode, "selected") == 0) {
        for (PetscInt i = 0; i < g_profiler_count; ++i) {
            if (g_profiler_registry[i].always_log && g_profiler_registry[i].total_call_count == 0) {
                LOG(GLOBAL, LOG_WARNING,
                    "profiling.timestep_output.functions lists '%s', which recorded no call in the first "
                    "step. It is not an instrumented function; the final summary lists the names that are.\n",
                    g_profiler_registry[i].name);
            }
        }
    }
 
    if (should_write && simCtx->rank == 0) {
        snprintf(filen, sizeof(filen), "%s/%s", simCtx->log_dir, simCtx->profilingTimestepFile);
        if (step == simCtx->StartStep + 1 && !simCtx->continueMode) {
            f = fopen(filen, "w");
            if (!f) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
            }
            PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
        } else {
            f = fopen(filen, "a");
            if (!f) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
            }
            if (step == simCtx->StartStep + 1 && ftell(f) == 0) {
                PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
            }
        }
        if (simCtx->continueMode && step == simCtx->StartStep + 1) {
            PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
        }
 
        for (PetscInt i = 0; i < g_profiler_count; ++i) {
            if (g_profiler_registry[i].current_step_call_count <= 0) {
                continue;
            }
            if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
                PetscFPrintf(
                    PETSC_COMM_SELF,
                    f,
                    "%d,%s,%lld,%.6f\n",
                    (int)step,
                    g_profiler_registry[i].name,
                    g_profiler_registry[i].current_step_call_count,
                    g_profiler_registry[i].current_step_time
                );
            }
        }
        fclose(f);
    }
 
    // Reset per-step counters for the next iteration
    for (PetscInt i = 0; i < g_profiler_count; ++i) {
        g_profiler_registry[i].current_step_time = 0.0;
        g_profiler_registry[i].current_step_call_count = 0;
    }
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref RuntimeMemoryLogSample().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see RuntimeMemoryLogSample()
 */
PetscErrorCode RuntimeMemoryLogSample(SimCtx *simCtx, PetscInt step, const char *event, const char *reason)
{
    PetscErrorCode ierr;
    PetscLogDouble process_current_bytes = 0.0;
    PetscLogDouble process_peak_bytes = 0.0;
    PetscLogDouble petsc_current_bytes = 0.0;
    PetscLogDouble petsc_peak_bytes = 0.0;
    PetscReal      local_values[5];
    PetscReal      global_values[5];
    PetscReal      process_current_mb = 0.0;
    PetscReal      process_peak_mb = 0.0;
    PetscReal      petsc_current_mb = 0.0;
    PetscReal      petsc_peak_mb = 0.0;
    PetscReal      process_change_mb = 0.0;
    char           path[(2 * PETSC_MAX_PATH_LEN) + 16];
    FILE           *f = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for RuntimeMemoryLogSample");
    if (!simCtx->runtimeMemoryLogEnabled) PetscFunctionReturn(0);
 
    ierr = PetscMemoryGetCurrentUsage(&process_current_bytes); CHKERRQ(ierr);
    ierr = PetscMemoryGetMaximumUsage(&process_peak_bytes); CHKERRQ(ierr);
    ierr = PetscMallocGetCurrentUsage(&petsc_current_bytes); CHKERRQ(ierr);
    ierr = PetscMallocGetMaximumUsage(&petsc_peak_bytes); CHKERRQ(ierr);
 
    process_current_mb = (PetscReal)(process_current_bytes / (1024.0 * 1024.0));
    process_peak_mb    = (PetscReal)(process_peak_bytes / (1024.0 * 1024.0));
    petsc_current_mb   = (PetscReal)(petsc_current_bytes / (1024.0 * 1024.0));
    petsc_peak_mb      = (PetscReal)(petsc_peak_bytes / (1024.0 * 1024.0));
    if (simCtx->runtimeMemoryLogHasPrevious) {
        process_change_mb = process_current_mb - simCtx->runtimeMemoryLogPreviousProcessMB;
    }
 
    local_values[0] = process_current_mb;
    local_values[1] = process_peak_mb;
    local_values[2] = petsc_current_mb;
    local_values[3] = petsc_peak_mb;
    local_values[4] = process_change_mb;
    ierr = MPI_Allreduce(local_values, global_values, 5, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    simCtx->runtimeMemoryLogPreviousProcessMB = process_current_mb;
    simCtx->runtimeMemoryLogHasPrevious = PETSC_TRUE;
 
    if (simCtx->rank == 0) {
        ierr = PetscSNPrintf(path, sizeof(path), "%s/%s", simCtx->log_dir, simCtx->runtimeMemoryLogFile); CHKERRQ(ierr);
        f = fopen(path, "a");
        if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open runtime memory log file: %s", path);
 
        if (!simCtx->runtimeMemoryLogStarted) {
            fprintf(f, "# PICurv runtime memory log\n");
            if (simCtx->continueMode) {
                fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
            }
            fprintf(
                f,
                "%-8s %-10s %22s %20s %22s %28s %22s %-18s\n",
                "Step",
                "Event",
                "Process Current MB Max",
                "Process Peak MB Max",
                "PETSc Allocated MB Max",
                "PETSc Peak Allocated MB Max",
                "Process Change MB Max",
                "Reason"
            );
            simCtx->runtimeMemoryLogStarted = PETSC_TRUE;
        }
 
        fprintf(
            f,
            "%-8" PetscInt_FMT " %-10s %22.3f %20.3f %22.3f %28.3f %22.3f %-18s\n",
            step,
            event ? event : "-",
            (double)global_values[0],
            (double)global_values[1],
            (double)global_values[2],
            (double)global_values[3],
            (double)global_values[4],
            (reason && reason[0]) ? reason : "-"
        );
        if ((event && (strcmp(event, "Shutdown") == 0 || strcmp(event, "Final") == 0))) {
            fflush(f);
        }
        fclose(f);
    }
 
    PetscFunctionReturn(0);
}
 
// Comparison function for qsort to sort by total_time in descending order
/**
 * @brief Order profiling records by their accumulated execution time.
 */
static int _CompareProfiledFunctions(const void *a, const void *b)
{
    const ProfiledFunction *func_a = (const ProfiledFunction *)a;
    const ProfiledFunction *func_b = (const ProfiledFunction *)b;
 
    if (func_a->total_time < func_b->total_time) return 1;
    if (func_a->total_time > func_b->total_time) return -1;
    return 0;
}
 
/**
 * @brief Implementation of \ref ProfilingFinalize().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ProfilingFinalize()
 */
PetscErrorCode ProfilingFinalize(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscInt rank = simCtx->rank;
    PetscFunctionBeginUser;
    if (!simCtx->profilingFinalSummary) PetscFunctionReturn(0);
    if (!rank) {
        
        char exec_mode_modifier[32] = "Unknown";
        if(simCtx->exec_mode == EXEC_MODE_SOLVER) PetscCall(PetscStrncpy(exec_mode_modifier, "Solver", sizeof(exec_mode_modifier)));
        else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) PetscCall(PetscStrncpy(exec_mode_modifier, "PostProcessor", sizeof(exec_mode_modifier)));
        //--- Step 0: Create a file viewer for log file
        FILE *f;
        char filen[PETSC_MAX_PATH_LEN + 128];
        ierr = PetscSNPrintf(filen, sizeof(filen), "%s/ProfilingSummary_%s.log",simCtx->log_dir,exec_mode_modifier); CHKERRQ(ierr);
 
        // Open the log file: append with section label in continue mode, truncate otherwise.
        if (simCtx->continueMode) {
            f = fopen(filen, "a");
            if (!f) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
            }
            fprintf(f, "\n=== Continuation from step %" PetscInt_FMT " ===\n", simCtx->StartStep);
        } else {
            f = fopen(filen, "w");
            if (!f) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
            }
        }
 
        // --- Step 1: Sort the data for readability ---
        qsort(g_profiler_registry, g_profiler_count, sizeof(ProfiledFunction), _CompareProfiledFunctions);
 
        // --- Step 2: Dynamically determine the width for the function name column ---
        PetscInt max_name_len = strlen("Function"); // Start with the header's length
        for (PetscInt i = 0; i < g_profiler_count; ++i) {
            if (g_profiler_registry[i].total_call_count > 0) {
                PetscInt len = strlen(g_profiler_registry[i].name);
                if (len > max_name_len) {
                    max_name_len = len;
                }
            }
        }
        // Add a little padding
        max_name_len += 2; 
 
        // --- Step 3: Define fixed widths for numeric columns for consistent alignment ---
        const int time_width = 18;
        const int count_width = 15;
        const int avg_width = 22;
 
        // --- Step 4: Print the formatted table ---
        PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
        PetscFPrintf(PETSC_COMM_SELF, f, "                         FINAL PROFILING SUMMARY (Sorted by Total Time)\n");
        PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
        
        // Header Row
        PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %-*s | %-*s | %-*s\n",
                    max_name_len, "Function",
                    time_width, "Total Time (s)",
                    count_width, "Call Count",
                    avg_width, "Avg. Time/Call (ms)");
        
        // Separator Line (dynamically sized)
        for (int i = 0; i < max_name_len; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
        PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
        for (int i = 0; i < time_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
        PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
        for (int i = 0; i < count_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
        PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
        for (int i = 0; i < avg_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
        PetscFPrintf(PETSC_COMM_SELF, f, "\n");
 
        // Data Rows
        for (PetscInt i = 0; i < g_profiler_count; ++i) {
            if (g_profiler_registry[i].total_call_count > 0) {
                double avg_time_ms = (g_profiler_registry[i].total_time / g_profiler_registry[i].total_call_count) * 1000.0;
                PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %*.*f | %*lld | %*.*f\n",
                            max_name_len, g_profiler_registry[i].name,
                            time_width, 6, g_profiler_registry[i].total_time,
                            count_width, g_profiler_registry[i].total_call_count,
                            avg_width, 6, avg_time_ms);
                PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------\n");
            }
        }
        PetscFPrintf(PETSC_COMM_SELF, f, "==================================================================================================================\n");
    
        fclose(f);
    }
 
    // --- Final Cleanup ---
    PetscFree(g_profiler_registry);
    g_profiler_registry = NULL;
    g_profiler_count = 0;
    g_profiler_capacity = 0;
    PetscFunctionReturn(0);
}
 
/*================================================================================*
 *                          PROGRESS BAR UTILITY                                  *
 *================================================================================*/
 
/**
 * @brief Internal helper implementation: `PrintProgressBar()`.
 * @details Local to this translation unit.
 */
void PrintProgressBar(PetscInt step, PetscInt startStep, PetscInt totalSteps, PetscReal currentTime)
{
    if (totalSteps <= 0) return;
 
    // --- Configuration ---
    const int barWidth = 50;
 
    // --- Calculation ---
    // Calculate progress as a fraction from 0.0 to 1.0
    PetscReal progress = (PetscReal)(step - startStep + 1) / totalSteps;
    // Ensure progress doesn't exceed 1.0 due to floating point inaccuracies
    if (progress > 1.0) progress = 1.0;
 
    int pos = (int)(barWidth * progress);
 
    // --- Printing ---
    // Carriage return moves cursor to the beginning of the line
    PetscPrintf(PETSC_COMM_SELF, "\rProgress: [");
 
    for (int i = 0; i < barWidth; ++i) {
        if (i < pos) {
            PetscPrintf(PETSC_COMM_SELF, "=");
        } else if (i == pos) {
            PetscPrintf(PETSC_COMM_SELF, ">");
        } else {
            PetscPrintf(PETSC_COMM_SELF, " ");
        }
    }
 
    // Print percentage, step count, and current time
    PetscPrintf(PETSC_COMM_SELF, "] %3d%% (Step %" PetscInt_FMT "/%" PetscInt_FMT ", t=%.4f)",
                (int)(progress * 100.0),
                step + 1,
                startStep + totalSteps,
                currentTime);
 
    // Flush the output buffer to ensure the bar is displayed immediately
    fflush(stdout);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_FIELD_MIN_MAX"
/**
 * @brief Implementation of \ref LOG_FIELD_MIN_MAX().
 * @details Full API contract is documented with the header declaration in `include/logging.h`.
 * @see LOG_FIELD_MIN_MAX()
 */
PetscErrorCode LOG_FIELD_MIN_MAX(UserCtx *user, FieldId field_id)
{
    PetscErrorCode ierr;
    PetscInt       i, j, k;
    DMDALocalInfo  info;
    
    FieldView      view;
    Vec            fieldVec = NULL;
    DM             dm = NULL;
    PetscInt       dof;
    FieldLayout    layout;
    const char    *fieldName = NULL;
    const char    *data_layout = NULL;
 
    PetscFunctionBeginUser;
 
    ierr = FieldGetView(user, field_id, &view); CHKERRQ(ierr);
    fieldName = view.descriptor->canonical_name;
    dm = view.dm;
    dof = view.descriptor->dof;
    layout = view.descriptor->layout;
    data_layout = FieldLayoutName(layout);
    fieldVec = (layout == FIELD_LAYOUT_COMPONENT_STAGGERED) ? view.local_vec : view.global_vec;
 
    ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
 
    // --- 2. Define Architecture-Aware Loop Bounds ---
    PetscInt i_start, i_end, j_start, j_end, k_start, k_end;
 
    if (layout == FIELD_LAYOUT_CELL_CENTERED) {
        // For cell-centered data, the physical values are stored from index 1 to N-1.
        // We find the intersection of the rank's owned range [xs, xe) with the
        // physical data range [1, IM-1).
        i_start = PetscMax(info.xs, 1);          i_end = PetscMin(info.xs + info.xm, user->IM);
        j_start = PetscMax(info.ys, 1);          j_end = PetscMin(info.ys + info.ym, user->JM);
        k_start = PetscMax(info.zs, 1);          k_end = PetscMin(info.zs + info.zm, user->KM);
    } else { // For Node- or Face-Centered data
        // The physical values are stored from index 0 to N-1.
        // We find the intersection of the rank's owned range [xs, xe) with the
        // physical data range [0, IM-1].
        i_start = PetscMax(info.xs, 0);          i_end = PetscMin(info.xs + info.xm, user->IM);
        j_start = PetscMax(info.ys, 0);          j_end = PetscMin(info.ys + info.ym, user->JM);
        k_start = PetscMax(info.zs, 0);          k_end = PetscMin(info.zs + info.zm, user->KM);
    }
 
    // --- 3. Barrier for clean, grouped output ---
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
    if (user->simCtx->rank == 0) {
        PetscPrintf(PETSC_COMM_SELF, "\n--- Field Ranges: [%s] (Layout: %s) ---\n", fieldName, data_layout);
    }
 
    // --- 4. Branch on DoF and perform calculation with correct bounds ---
    if (dof == 1) {
        PetscReal    localMin = PETSC_MAX_REAL, localMax = PETSC_MIN_REAL;
        PetscReal    globalMin, globalMax;
        const PetscScalar  ***array;
 
        ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
        for (k = k_start; k < k_end; k++) {
            for (j = j_start; j < j_end; j++) {
                for (i = i_start; i < i_end; i++) {
                    localMin = PetscMin(localMin, array[k][j][i]);
                    localMax = PetscMax(localMax, array[k][j][i]);
                }
            }
        }
        ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
 
        ierr = MPI_Allreduce(&localMin, &globalMin, 1, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(&localMax, &globalMax, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  [Rank %d] Local Range:  [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin, localMax);
        ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
        if (user->simCtx->rank == 0) {
           PetscPrintf(PETSC_COMM_SELF, "  Global Range:         [ %11.4e , %11.4e ]\n", globalMin, globalMax);
        }
 
    } else if (dof == 3) {
        Cmpnts localMin = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
        Cmpnts localMax = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
        Cmpnts globalMin, globalMax;
        const Cmpnts ***array;
 
        ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
        for (k = k_start; k < k_end; k++) {
            for (j = j_start; j < j_end; j++) {
                for (i = i_start; i < i_end; i++) {
                    localMin.x = PetscMin(localMin.x, array[k][j][i].x);
                    localMin.y = PetscMin(localMin.y, array[k][j][i].y);
                    localMin.z = PetscMin(localMin.z, array[k][j][i].z);
                    localMax.x = PetscMax(localMax.x, array[k][j][i].x);
                    localMax.y = PetscMax(localMax.y, array[k][j][i].y);
                    localMax.z = PetscMax(localMax.z, array[k][j][i].z);
                }
            }
        }
        ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
 
        ierr = MPI_Allreduce(&localMin, &globalMin, 3, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(&localMax, &globalMax, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  [Rank %d] Local X-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.x, localMax.x);
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  [Rank %d] Local Y-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.y, localMax.y);
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  [Rank %d] Local Z-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.z, localMax.z);
        ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
 
        if (user->simCtx->rank == 0) {
            PetscPrintf(PETSC_COMM_SELF, "  [Global] X-Range: [ %11.4e , %11.4e ]\n", globalMin.x, globalMax.x);
            PetscPrintf(PETSC_COMM_SELF, "  [Global] Y-Range: [ %11.4e , %11.4e ]\n", globalMin.y, globalMax.y);
            PetscPrintf(PETSC_COMM_SELF, "  [Global] Z-Range: [ %11.4e , %11.4e ]\n", globalMin.z, globalMax.z);
        }
 
    } else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "LogFieldStatistics only supports fields with 1 or 3 components, but field '%s' has %" PetscInt_FMT ".", fieldName, dof);
    }
 
    // --- 5. Final barrier for clean output ordering ---
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
    if (user->simCtx->rank == 0) {
        PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------\n\n");
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LogFieldAnatomyView"
/**
 * @brief Shared architecture-aware anatomy logger for catalog and transient fields.
 */
static PetscErrorCode LogFieldAnatomyView(UserCtx *user, const char *field_name,
                                          const char *stage_name, DM dm,
                                          Vec vec_local, PetscInt dof,
                                          FieldLayout layout)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    const char    *data_layout = FieldLayoutName(layout);
    char           dominant_dir = '\0'; // 'x', 'y', 'z' for face-centered, 'm' for mixed (Ucont)
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    PetscCheck(user != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
    PetscCheck(field_name != NULL && stage_name != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Field and stage labels cannot be NULL.");
    PetscCheck(dm != NULL && vec_local != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "Field '%s' has no available DM/local vector for anatomy logging.", field_name);
    PetscCheck(dof == 1 || dof == 3, PETSC_COMM_SELF, PETSC_ERR_SUP,
               "Field anatomy logging supports one- or three-component fields; '%s' has %d.",
               field_name, dof);
 
    if (layout == FIELD_LAYOUT_I_FACE) dominant_dir = 'x';
    else if (layout == FIELD_LAYOUT_J_FACE) dominant_dir = 'y';
    else if (layout == FIELD_LAYOUT_K_FACE) dominant_dir = 'z';
    else if (layout == FIELD_LAYOUT_COMPONENT_STAGGERED) dominant_dir = 'm';
 
    // --- 2. Get Grid Info and Array Pointers ---
    ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
 
    ierr = PetscBarrier(NULL);
    PetscPrintf(PETSC_COMM_WORLD, "\n--- Field Anatomy Log: [%s] | Stage: [%s] | Layout: [%s] ---\n", field_name, stage_name, data_layout);
 
    // Global physical dimensions (number of cells)
    PetscInt im_phys = user->IM;
    PetscInt jm_phys = user->JM;
    PetscInt km_phys = user->KM;
 
    // Slice through the center of the local domain
    PetscInt i_mid = (PetscInt)(info.xs + 0.5 * info.xm) - 1;
    PetscInt j_mid = (PetscInt)(info.ys + 0.5 * info.ym) - 1;
    PetscInt k_mid = (PetscInt)(info.zs + 0.5 * info.zm) - 1;
 
    // --- 3. Print Boundary Information based on Data Layout ---
 
    // ======================================================================
    // === CASE 1: Cell-Centered Fields (Ucat, P) - USES SHIFTED INDEX    ===
    // ======================================================================
    if (layout == FIELD_LAYOUT_CELL_CENTERED) {
        const void *l_arr;
        ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
 
 
        // --- I-Direction Boundaries ---
        if (info.xs == 0) { // Rank on -Xi boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][0])      = ", rank, 0);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][0]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][0].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][0].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][0].z);
 
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][0])      = ", rank, 1);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][1]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][1].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][1].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][1].z);
        }
        if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][%d]) = ", rank, im_phys - 1, im_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys - 1]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].z);
 
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][%d]) = ", rank, im_phys, im_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].z);
        } 
 
        // --- J-Direction Boundaries ---
        if (info.ys == 0) { // Rank on -Eta boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][0][i])      = ", rank, 0);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][0][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][0][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][0][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][0][i_mid].z);
            
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][0][i])      = ", rank, 1);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][1][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][1][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][1][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][1][i_mid].z);
        }
 
        if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][%d][i]) = ", rank, jm_phys - 1, jm_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys - 1][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].z);
 
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][%d][i]) = ", rank, jm_phys, jm_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].z);
        }
 
        // --- K-Direction Boundaries ---
        if (info.zs == 0) { // Rank on -Zeta boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[0][j][i])      = ", rank, 0);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[0][j_mid][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[0][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[0][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[0][j_mid][i_mid].z);
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[0][j][i])      = ", rank, 1);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[1][j_mid][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[1][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[1][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[1][j_mid][i_mid].z);
        }
        if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[%d][j][i]) = ", rank, km_phys - 1, km_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys - 1][j_mid][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].z);
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[%d][j][i]) = ", rank, km_phys, km_phys - 2);
            if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys][j_mid][i_mid]);
            else       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].z);
        }
        ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
    }
    // ======================================================================
    // === CASE 2: Face-Centered Fields - NUANCED DIRECTIONAL LOGIC       ===
    // ======================================================================
    else if (layout == FIELD_LAYOUT_I_FACE ||
             layout == FIELD_LAYOUT_J_FACE ||
             layout == FIELD_LAYOUT_K_FACE ||
             layout == FIELD_LAYOUT_COMPONENT_STAGGERED) {
        const Cmpnts ***l_arr;
        ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
 
        // --- I-Direction Boundaries ---
        if (info.xs == 0) { // Rank on -Xi boundary
            if (dominant_dir == 'x') { // Node-like in I-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (First Phys. X-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
            } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][0]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][0]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[k_mid][j_mid][1].x, l_arr[k_mid][j_mid][1].y, l_arr[k_mid][j_mid][1].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (1st X-Face) = %.5f\n", rank, 0, l_arr[k_mid][j_mid][0].x);
            }
        }
        if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
            if (dominant_dir == 'x') { // Node-like in I-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Last Phys. X-Face)  = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys-1].y, l_arr[k_mid][j_mid][im_phys - 1].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost Location)     = (%.5f, %.5f, %.5f)\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
            } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][%d]) = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, im_phys - 2, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys - 1].y, l_arr[k_mid][j_mid][im_phys-1].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][%d]) = (%.5f, %.5f, %.5f)\n", rank, im_phys, im_phys - 2, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (Last X-Face)   = %.5f\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (Ghost Location)    = %.5f\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x);
            }
        }
 
        // --- J-Direction Boundaries ---
        if (info.ys == 0) { // Rank on -Eta boundary
            if (dominant_dir == 'y') { // Node-like in J-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (First Phys. Y-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
            } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][0][i]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][0][i]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[k_mid][1][i_mid].x, l_arr[k_mid][1][i_mid].y, l_arr[k_mid][1][i_mid].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (1st Y-Face) = %.5f\n", rank, 0, l_arr[k_mid][0][i_mid].y);
            }
        }
        if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
             if (dominant_dir == 'y') { // Node-like in J-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Last Phys. Y-Face)  = (%.5f, %.5f, %.5f)\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost Location)     = (%.5f, %.5f, %.5f)\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
            } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][%d][i]) = (%.5f, %.5f, %.5f)\n", rank, jm_phys-1, jm_phys-2, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][%d][i]) = (%.5f, %.5f, %.5f)\n", rank, jm_phys, jm_phys-2, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (Last Y-Face)    = %.5f\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].y); 
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (Ghost Location)   = %.5f\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].y);
            }
        }
 
        // --- K-Direction Boundaries ---
        if (info.zs == 0) { // Rank on -Zeta boundary
            if (dominant_dir == 'z') { // Node-like in K-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (First Phys. Z-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
            } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[0][j][i]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[0][j][i]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[1][j_mid][i_mid].x, l_arr[1][j_mid][i_mid].y, l_arr[1][j_mid][i_mid].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (1st Z-Face) = %.5f\n", rank, 0, l_arr[0][j_mid][i_mid].z);
            }
        }
        if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
             if (dominant_dir == 'z') { // Node-like in K-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Last Phys. Z-Face)  = (%.5f, %.5f, %.5f)\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].x, l_arr[km_phys - 1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Ghost Location)     = (%.5f, %.5f, %.5f)\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
            } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Value for Cell[%d][j][i]) = (%.5f, %.5f, %.5f)\n", rank, km_phys-1, km_phys-2, l_arr[km_phys-1][j_mid][i_mid].x, l_arr[km_phys-1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Ghost for Cell[%d][j][i]) = (%.5f, %.5f, %.5f)\n", rank, km_phys, km_phys-2, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
            } else if (dominant_dir == 'm') { // Ucont: Mixed
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (Last Z-Face) = %.5f\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].z);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (Ghost Loc.) = %.5f\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].z);
 
            }
        }
        ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
    }
    // ======================================================================
    // === CASE 3: Node-Centered Fields - USES DIRECT INDEX               ===
    // ======================================================================
    else if (layout == FIELD_LAYOUT_NODE_CENTERED) {
        if (dof == 1) {
            const PetscReal ***l_arr;
            ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
 
            if (info.xs == 0)
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (First Phys. Node) = %.5f\n", rank, 0, l_arr[k_mid][j_mid][0]);
            if (info.xs + info.xm == info.mx) {
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Last Phys. Node)  = %.5f\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1]);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Unused/Ghost Loc) = %.5f\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys]);
            }
            if (info.ys == 0)
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (First Phys. Node) = %.5f\n", rank, 0, l_arr[k_mid][0][i_mid]);
            if (info.ys + info.ym == info.my) {
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Last Phys. Node)  = %.5f\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid]);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Unused/Ghost Loc) = %.5f\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid]);
            }
            if (info.zs == 0)
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (First Phys. Node) = %.5f\n", rank, 0, l_arr[0][j_mid][i_mid]);
            if (info.zs + info.zm == info.mz) {
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Last Phys. Node)  = %.5f\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid]);
                PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Unused/Ghost Loc) = %.5f\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid]);
            }
            ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
        } else {
            const Cmpnts ***l_arr;
            ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
 
        // --- I-Direction Boundaries ---
        if (info.xs == 0) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
        }
        if (info.xs + info.xm == info.mx) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Last Phys. Node)  = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys - 1].y, l_arr[k_mid][j_mid][im_phys - 1].z);
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
        }
        // --- J-Direction Boundaries ---
        if (info.ys == 0) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
        }
        if (info.ys + info.ym == info.my) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Last Phys. Node)  = (%.5f, %.5f, %.5f)\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
        }
        // --- K-Direction Boundaries ---
        if (info.zs == 0) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
        }
        if(info.zs + info.zm == info.mz) {
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Last Phys. Node)  = (%.5f, %.5f, %.5f)\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].x, l_arr[km_phys - 1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
            PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
        }
            ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
        }
    }
    else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "LOG_FIELD_ANATOMY encountered unsupported layout %d for field '%s'.",
                (int)layout, field_name);
    }
 
    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
    ierr = PetscBarrier(NULL);
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_FIELD_ANATOMY"
/**
 * @brief Resolves a persistent field view and emits its layout-aware boundary anatomy.
 * @see LOG_FIELD_ANATOMY()
 */
PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, FieldId field_id, const char *stage_name)
{
    FieldView view;
 
    PetscFunctionBeginUser;
    PetscCall(FieldGetView(user, field_id, &view));
    PetscCall(LogFieldAnatomyView(user, view.descriptor->canonical_name, stage_name,
                                  view.dm, view.local_vec, view.descriptor->dof,
                                  view.descriptor->layout));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref IsStatisticsConsoleSnapshotEnabled().
 * @see IsStatisticsConsoleSnapshotEnabled()
 */
PetscBool IsStatisticsConsoleSnapshotEnabled(const SimCtx *simCtx)
{
    if (!FieldStatisticsIsActive(simCtx)) return PETSC_FALSE;
    return (PetscBool)(simCtx->statisticsConsoleOutputFreq > 0 && get_log_level() >= LOG_INFO);
}
 
/**
 * @brief Implementation of \ref ShouldEmitPeriodicStatisticsConsoleSnapshot().
 * @see ShouldEmitPeriodicStatisticsConsoleSnapshot()
 */
PetscBool ShouldEmitPeriodicStatisticsConsoleSnapshot(const SimCtx *simCtx, PetscInt completed_step)
{
    if (!IsStatisticsConsoleSnapshotEnabled(simCtx)) return PETSC_FALSE;
    if (completed_step < 0) return PETSC_FALSE;
    return (PetscBool)((completed_step % simCtx->statisticsConsoleOutputFreq) == 0);
}
 
/**
 * @brief Implementation of \ref EmitStatisticsConsoleSnapshot().
 * @see EmitStatisticsConsoleSnapshot()
 */
PetscErrorCode EmitStatisticsConsoleSnapshot(UserCtx *user, const SimCtx *simCtx, PetscInt step)
{
    PetscFunctionBeginUser;
    PetscCheck(simCtx != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL.");
    if (!IsStatisticsConsoleSnapshotEnabled(simCtx)) PetscFunctionReturn(0);
 
    LOG(GLOBAL, LOG_INFO, "Statistics windows at step %d (%d window(s)):\n",
        step, simCtx->fieldStatisticsWindowCount);
    for (PetscInt w = 0; w < simCtx->fieldStatisticsWindowCount; ++w) {
        const PicurvWindow *window = &simCtx->fieldStatisticsWindows[w];
        char progress[16];
        char coverage[32] = "";
 
        if (window->definition.bounded) {
            PetscCall(PetscSNPrintf(progress, sizeof(progress), "%5.1f%%",
                                    100.0 * (double)PicurvWindowProgress(window)));
        } else {
            PetscCall(PetscSNPrintf(progress, sizeof(progress), "  open"));
        }
 
        /* Mask health: with a moving body different points see different numbers of
         * states, and the spread of per-point valid fraction is what tells an
         * operator whether the window covers the domain evenly. Reported only on the
         * first block, since the reduction it performs is already collective. */
        if (user && user->fieldStatisticsStorage && window->sample_count > 0) {
            PetscReal lowest = 1.0, highest = 0.0;
 
            PetscCall(PicurvWindowValidFractionRange(user, &window->definition,
                                                     &user->fieldStatisticsStorage[w],
                                                     window->sample_count, &lowest, &highest));
            PetscCall(PetscSNPrintf(coverage, sizeof(coverage), " valid=[%.2f,%.2f]",
                                    (double)lowest, (double)highest));
        }
 
        LOG(GLOBAL, LOG_INFO,
            "  %-24s %-8s samples=%-8d weight=%-12.6g represented=%-12.6g progress=%s%s\n",
            window->definition.name, PicurvWindowStateName(window->state),
            window->sample_count, (double)window->total_weight,
            (double)window->represented_time, progress, coverage);
    }
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_CORNER_FIELD_ANATOMY"
/**
 * @brief Emits anatomy for the transient scalar or vector corner-staging field.
 * @see LOG_CORNER_FIELD_ANATOMY()
 */
PetscErrorCode LOG_CORNER_FIELD_ANATOMY(UserCtx *user, FieldId corner_field_id, const char *stage_name)
{
    FieldView view;
 
    PetscFunctionBeginUser;
    PetscCheck(user != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
    PetscCheck(corner_field_id == FIELD_ID_CELL_SCALAR_AT_CORNER ||
               corner_field_id == FIELD_ID_CELL_VECTOR_AT_CORNER,
               PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Corner anatomy logging expects a corner-staging field identity.");
    /* The caller states which workspace it used, so the degree of freedom and DM
     * come from the catalog instead of being inferred from a cached vector. */
    PetscCall(FieldGetView(user, corner_field_id, &view));
    PetscCall(LogFieldAnatomyView(user, view.descriptor->canonical_name, stage_name, view.dm,
                                  view.local_vec, view.descriptor->dof,
                                  view.descriptor->layout));
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_INTERPOLATION_ERROR"
/**
 * @brief Implementation of \ref LOG_INTERPOLATION_ERROR().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_INTERPOLATION_ERROR()
 */
PetscErrorCode LOG_INTERPOLATION_ERROR(UserCtx *user)
{
    SimCtx *simCtx = user->simCtx;
    PetscErrorCode ierr;
    DM swarm = user->swarm;
    Vec positionVec, analyticalvelocityVec, velocityVec, errorVec;
    PetscReal Interpolation_error = 0.0;
    PetscReal Maximum_Interpolation_error = 0.0;
    PetscReal AnalyticalSolution_magnitude = 0.0;
    PetscReal ErrorPercentage = 0.0;
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating global vectors.\n");
    ierr = DMSwarmCreateGlobalVectorFromField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), &positionVec); CHKERRQ(ierr);
    ierr = DMSwarmCreateGlobalVectorFromField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), &velocityVec); CHKERRQ(ierr);
    
    ierr = VecDuplicate(positionVec, &analyticalvelocityVec); CHKERRQ(ierr);
    ierr = VecCopy(positionVec, analyticalvelocityVec); CHKERRQ(ierr);
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing analytical solution.\n");
    ierr = SetAnalyticalSolutionForParticles(analyticalvelocityVec, simCtx); CHKERRQ(ierr);
    
    ierr = VecDuplicate(analyticalvelocityVec, &errorVec); CHKERRQ(ierr);
    ierr = VecCopy(analyticalvelocityVec, errorVec); CHKERRQ(ierr);
    
    ierr = VecNorm(analyticalvelocityVec, NORM_2, &AnalyticalSolution_magnitude); CHKERRQ(ierr);
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing error.\n");
    ierr = VecAXPY(errorVec, -1.0, velocityVec); CHKERRQ(ierr);
    ierr = VecNorm(errorVec, NORM_2, &Interpolation_error); CHKERRQ(ierr);
    ierr = VecNorm(errorVec,NORM_INFINITY,&Maximum_Interpolation_error); CHKERRQ(ierr);
    
    ErrorPercentage = (AnalyticalSolution_magnitude > 0) ?
                      (Interpolation_error / AnalyticalSolution_magnitude * 100.0) : 0.0;
 
    /* --- CSV output (always, rank 0 only) --- */
    if (simCtx->rank == 0) {
        char csv_path[PETSC_MAX_PATH_LEN + 32];
        ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/interpolation_error.csv", simCtx->analysis_dir); CHKERRQ(ierr);
        FILE *f = fopen(csv_path, "a");
        if (f) {
            if (ftell(f) == 0) {
                fprintf(f, "step,time,L2_error,Linf_error,L2_analytical,error_pct,physical_time\n");
            }
            if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
                fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
            }
            /* ti is solver time, not a step count; physical_time is it in seconds. */
            PetscReal t = simCtx->ti;
            PetscReal physical_time = 0.0;
            ierr = PicurvPhysicalTime(simCtx, t, &physical_time); CHKERRQ(ierr);
            fprintf(f, "%d,%.6e,%.6e,%.6e,%.6e,%.4f,%.6e\n",
                    (int)simCtx->step, t,
                    Interpolation_error, Maximum_Interpolation_error,
                    AnalyticalSolution_magnitude, ErrorPercentage, physical_time);
            fclose(f);
        }
    }
 
    /* --- Console output (only at INFO level or above) --- */
    if (get_log_level() >= LOG_INFO) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Interpolation error (%%): %g\n", ErrorPercentage);
        PetscPrintf(PETSC_COMM_WORLD, "Interpolation error (%%): %g\n", ErrorPercentage);
        LOG_ALLOW(GLOBAL, LOG_INFO, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
        PetscPrintf(PETSC_COMM_WORLD, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
    }
 
    ierr = VecDestroy(&analyticalvelocityVec); CHKERRQ(ierr);
    ierr = VecDestroy(&errorVec); CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), &positionVec); CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), &velocityVec); CHKERRQ(ierr);
    
    return 0;
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_SCATTER_METRICS"
/**
 * @brief Implementation of \ref LOG_SCATTER_METRICS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_SCATTER_METRICS()
 */
PetscErrorCode LOG_SCATTER_METRICS(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx        *simCtx = NULL;
    DMDALocalInfo  info;
    PetscInt       xs, xe, ys, ye, zs, ze, mx, my, mz;
    PetscInt       lxs, lxe, lys, lye, lzs, lze;
    Vec            reference_vec = NULL;
    PetscReal    ***psi = NULL;
    PetscReal    ***psi_ref = NULL;
    PetscReal    ***aj = NULL;
    PetscReal    ***count = NULL;
    PetscReal     *particle_psi = NULL;
    PetscInt       nlocal = 0;
    PetscReal      local_l1 = 0.0, global_l1 = 0.0;
    PetscReal      local_l2_sq = 0.0, global_l2_sq = 0.0;
    PetscReal      local_linf = 0.0, global_linf = 0.0;
    PetscReal      local_ref_l2_sq = 0.0, global_ref_l2_sq = 0.0;
    PetscReal      local_grid_integral = 0.0, global_grid_integral = 0.0;
    PetscReal      local_domain_volume = 0.0, global_domain_volume = 0.0;
    PetscReal      local_particle_sum = 0.0, global_particle_sum = 0.0;
    PetscInt64     local_particle_count = 0, global_particle_count = 0;
    PetscInt64     local_cell_count = 0, global_cell_count = 0;
    PetscInt64     local_occupied_count = 0, global_occupied_count = 0;
    PetscReal      particle_integral = 0.0;
    PetscReal      occupancy_fraction = 0.0;
    PetscReal      mean_particles_per_occupied_cell = 0.0;
    PetscReal      l2_error = 0.0;
    PetscReal      relative_l2_error = 0.0;
 
    PetscFunctionBeginUser;
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
    simCtx = user->simCtx;
    if (!VerificationScalarOverrideActive(simCtx) || !user->swarm || !user->Psi || !user->ParticleCount) {
        PetscFunctionReturn(0);
    }
 
    info = user->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;
    mx = info.mx; my = info.my; mz = info.mz;
    lxs = (xs == 0) ? xs + 1 : xs; lxe = (xe == mx) ? xe - 1 : xe;
    lys = (ys == 0) ? ys + 1 : ys; lye = (ye == my) ? ye - 1 : ye;
    lzs = (zs == 0) ? zs + 1 : zs; lze = (ze == mz) ? ze - 1 : ze;
 
    ierr = VecDuplicate(user->Psi, &reference_vec); CHKERRQ(ierr);
    ierr = SetAnalyticalScalarFieldAtCellCenters(user, reference_vec); CHKERRQ(ierr);
 
    ierr = DMDAVecGetArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
 
    for (PetscInt k = lzs; k < lze; ++k) {
        for (PetscInt j = lys; j < lye; ++j) {
            for (PetscInt i = lxs; i < lxe; ++i) {
                const PetscReal cell_volume = (PetscAbsReal(aj[k][j][i]) > 1.0e-14) ? (1.0 / aj[k][j][i]) : 0.0;
                const PetscReal err = psi[k][j][i] - psi_ref[k][j][i];
                local_cell_count += 1;
                local_domain_volume += cell_volume;
                local_grid_integral += psi[k][j][i] * cell_volume;
                local_l1 += PetscAbsReal(err) * cell_volume;
                local_l2_sq += err * err * cell_volume;
                local_ref_l2_sq += psi_ref[k][j][i] * psi_ref[k][j][i] * cell_volume;
                local_linf = PetscMax(local_linf, PetscAbsReal(err));
                if (count[k][j][i] > 0.0) local_occupied_count += 1;
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
    ierr = VecDestroy(&reference_vec); CHKERRQ(ierr);
 
    ierr = DMSwarmGetLocalSize(user->swarm, &nlocal); CHKERRQ(ierr);
    local_particle_count = (PetscInt64)nlocal;
    if (nlocal > 0) {
        ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI), NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
        for (PetscInt p = 0; p < nlocal; ++p) local_particle_sum += particle_psi[p];
        ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI), NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
    }
 
    ierr = MPI_Allreduce(&local_l1, &global_l1, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_l2_sq, &global_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_linf, &global_linf, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_ref_l2_sq, &global_ref_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_grid_integral, &global_grid_integral, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_domain_volume, &global_domain_volume, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_particle_sum, &global_particle_sum, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_particle_count, &global_particle_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_cell_count, &global_cell_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&local_occupied_count, &global_occupied_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    l2_error = PetscSqrtReal(global_l2_sq);
    relative_l2_error = (global_ref_l2_sq > 0.0) ? (l2_error / PetscSqrtReal(global_ref_l2_sq)) : 0.0;
    occupancy_fraction = (global_cell_count > 0) ? ((PetscReal)global_occupied_count / (PetscReal)global_cell_count) : 0.0;
    mean_particles_per_occupied_cell =
      (global_occupied_count > 0) ? ((PetscReal)global_particle_count / (PetscReal)global_occupied_count) : 0.0;
    particle_integral =
      (global_particle_count > 0) ? (global_domain_volume * global_particle_sum / (PetscReal)global_particle_count) : 0.0;
 
    if (simCtx->rank == 0) {
        char csv_path[PETSC_MAX_PATH_LEN + 32];
        FILE *f = NULL;
        ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/scatter_metrics.csv", simCtx->analysis_dir); CHKERRQ(ierr);
        f = fopen(csv_path, "a");
        if (f) {
            if (ftell(f) == 0) {
                fprintf(f,
                        "step,time,total_particles,total_cells,occupied_cells,occupancy_fraction,"
                        "mean_particles_per_occupied_cell,particle_integral,grid_integral,"
                        "conservation_error_abs,L1_error,L2_error,Linf_error,relative_L2_error,"
                        "physical_time\n");
            }
            if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
                fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
            }
            PetscReal physical_time = 0.0;
            ierr = PicurvPhysicalTime(simCtx, simCtx->ti, &physical_time); CHKERRQ(ierr);
            fprintf(f, "%d,%.6e,%lld,%lld,%lld,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e\n",
                    (int)simCtx->step,
                    (double)simCtx->ti,
                    (long long)global_particle_count,
                    (long long)global_cell_count,
                    (long long)global_occupied_count,
                    (double)occupancy_fraction,
                    (double)mean_particles_per_occupied_cell,
                    (double)particle_integral,
                    (double)global_grid_integral,
                    (double)PetscAbsReal(global_grid_integral - particle_integral),
                    (double)global_l1,
                    (double)l2_error,
                    (double)global_linf,
                    (double)relative_l2_error,
                    (double)physical_time);
            fclose(f);
        }
    }
 
    if (get_log_level() >= LOG_INFO) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter relative L2 error: %.6e\n", (double)relative_l2_error);
        LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter occupancy fraction: %.6e\n", (double)occupancy_fraction);
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ResetSearchMetrics"
/**
 * @brief Implementation of \ref ResetSearchMetrics().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see ResetSearchMetrics()
 */
PetscErrorCode ResetSearchMetrics(SimCtx *simCtx)
{
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL for ResetSearchMetrics.");
 
    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;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_SEARCH_METRICS"
/**
 * @brief Implementation of \ref LOG_SEARCH_METRICS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_SEARCH_METRICS()
 */
PetscErrorCode LOG_SEARCH_METRICS(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx         *simCtx = NULL;
    PetscInt       totalParticles = 0;
    PetscReal      local_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
    PetscReal      global_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
    PetscReal      meanTraversalSteps = 0.0;
    PetscReal      searchFailureFraction = 0.0;
    PetscReal      searchWorkIndex = 0.0;
    PetscReal      reSearchFraction = 0.0;
    long long      searchAttempts = 0;
    long long      searchPopulation = 0;
    long long      searchLocatedCount = 0;
    long long      searchLostCount = 0;
    long long      traversalStepsSum = 0;
    long long      reSearchCount = 0;
    long long      tieBreakCount = 0;
    long long      boundaryClampCount = 0;
    long long      bboxGuessSuccessCount = 0;
    long long      bboxGuessFallbackCount = 0;
    long long      maxTraversalFailCount = 0;
    long long      maxTraversalSteps = 0;
    long long      maxPassDepth = 0;
    MPI_Op         reduction_op = MPI_OP_NULL;
 
    PetscFunctionBeginUser;
    if (!user || !user->simCtx) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx and SimCtx are required for LOG_SEARCH_METRICS.");
    }
    simCtx = user->simCtx;
 
    if (simCtx->np <= 0) {
        PetscFunctionReturn(0);
    }
 
    ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
 
    local_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] = (PetscReal)simCtx->searchMetrics.searchAttempts;
    local_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] = (PetscReal)simCtx->searchMetrics.searchPopulation;
    local_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] = (PetscReal)simCtx->searchMetrics.searchLocatedCount;
    local_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] = (PetscReal)simCtx->searchMetrics.searchLostCount;
    local_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.traversalStepsSum;
    local_metrics[SEARCH_METRIC_SUM_RESEARCH] = (PetscReal)simCtx->searchMetrics.reSearchCount;
    local_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] = (PetscReal)simCtx->searchMetrics.tieBreakCount;
    local_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] = (PetscReal)simCtx->searchMetrics.boundaryClampCount;
    local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] = (PetscReal)simCtx->searchMetrics.bboxGuessSuccessCount;
    local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] = (PetscReal)simCtx->searchMetrics.bboxGuessFallbackCount;
    local_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] = (PetscReal)simCtx->searchMetrics.maxTraversalFailCount;
    local_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.maxTraversalSteps;
    local_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] = (PetscReal)simCtx->searchMetrics.maxParticlePassDepth;
 
    ierr = MPI_Op_create(SearchMetricsReduceOp, PETSC_TRUE, &reduction_op); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(local_metrics, global_metrics, SEARCH_METRIC_REDUCTION_LEN, MPIU_REAL, reduction_op, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Op_free(&reduction_op); CHKERRMPI(ierr);
    reduction_op = MPI_OP_NULL;
 
    searchAttempts = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] + 0.5);
    searchPopulation = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] + 0.5);
    searchLocatedCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] + 0.5);
    searchLostCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] + 0.5);
    traversalStepsSum = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] + 0.5);
    reSearchCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_RESEARCH] + 0.5);
    tieBreakCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] + 0.5);
    boundaryClampCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] + 0.5);
    bboxGuessSuccessCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] + 0.5);
    bboxGuessFallbackCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] + 0.5);
    maxTraversalFailCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] + 0.5);
    maxTraversalSteps = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] + 0.5);
    maxPassDepth = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] + 0.5);
 
    if (searchAttempts > 0) {
        meanTraversalSteps = (PetscReal)traversalStepsSum / (PetscReal)searchAttempts;
    }
    if (searchPopulation > 0) {
        searchFailureFraction = (PetscReal)searchLostCount / (PetscReal)searchPopulation;
        searchWorkIndex = (PetscReal)traversalStepsSum / (PetscReal)searchPopulation;
        reSearchFraction = (PetscReal)reSearchCount / (PetscReal)searchPopulation;
    }
 
    if (simCtx->rank == 0) {
        char csv_path[PETSC_MAX_PATH_LEN + 32];
        FILE *f = NULL;
        PetscReal searchPhysicalTime = 0.0;
 
        ierr = PicurvPhysicalTime(simCtx, simCtx->ti, &searchPhysicalTime); CHKERRQ(ierr);
        ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/search_metrics.csv", simCtx->analysis_dir); CHKERRQ(ierr);
        f = fopen(csv_path, "a");
        if (!f) {
            LOG_ALLOW(GLOBAL, LOG_WARNING, "LOG_SEARCH_METRICS: could not open '%s' for writing.\n", csv_path);
        } else {
            if (ftell(f) == 0) {
                fprintf(f,
                        "step,time,total_particles,lost,lost_cumulative,migrated,migration_passes,search_attempts,"
                        "mean_traversal_steps,max_traversal_steps,tie_break_count,boundary_clamp_count,"
                        "bbox_guess_success_count,bbox_guess_fallback_count,max_particle_pass_depth,load_imbalance,"
                        "search_population,search_located_count,search_lost_count,traversal_steps_sum,re_search_count,"
                        "max_traversal_fail_count,search_failure_fraction,search_work_index,re_search_fraction,"
                        "physical_time,lost_psi_sum\n");
            }
            if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
                fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
            }
            fprintf(f,
                    "%d,%.6e,%d,%d,%d,%d,%d,%lld,%.6e,%lld,%lld,%lld,%lld,%lld,%lld,%.6e,%lld,%lld,%lld,%lld,%lld,%lld,%.6e,%.6e,%.6e,%.6e,%.10e\n",
                    (int)simCtx->step,
                    (double)simCtx->ti,
                    (int)totalParticles,
                    (int)simCtx->particlesLostLastStep,
                    (int)simCtx->particlesLostCumulative,
                    (int)simCtx->particlesMigratedLastStep,
                    (int)simCtx->migrationPassesLastStep,
                    searchAttempts,
                    (double)meanTraversalSteps,
                    maxTraversalSteps,
                    tieBreakCount,
                    boundaryClampCount,
                    bboxGuessSuccessCount,
                    bboxGuessFallbackCount,
                    maxPassDepth,
                    (double)simCtx->particleLoadImbalance,
                    searchPopulation,
                    searchLocatedCount,
                    searchLostCount,
                    traversalStepsSum,
                    reSearchCount,
                    maxTraversalFailCount,
                    (double)searchFailureFraction,
                    (double)searchWorkIndex,
                    (double)reSearchFraction,
                    (double)searchPhysicalTime,
                    (double)simCtx->particlesLostScalarLastStep);
            fclose(f);
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Search metrics: sff=%.3e swi=%.3e re_search=%.3e lost(step/total)=%d/%d migrated=%d passes=%d traversal(mean/max)=%.2f/%lld tie_breaks=%lld max_pass_depth=%lld\n",
              (double)searchFailureFraction,
              (double)searchWorkIndex,
              (double)reSearchFraction,
              (int)simCtx->particlesLostLastStep,
              (int)simCtx->particlesLostCumulative,
              (int)simCtx->particlesMigratedLastStep,
              (int)simCtx->migrationPassesLastStep,
              (double)meanTraversalSteps,
              maxTraversalSteps,
              tieBreakCount,
              maxPassDepth);
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "CalculateAdvancedParticleMetrics"
/**
 * @brief Internal helper implementation: `CalculateAdvancedParticleMetrics()`.
 * @details Local to this translation unit.
 */
PetscErrorCode CalculateAdvancedParticleMetrics(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
    PetscMPIInt    size, rank;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- 1. Particle Load Imbalance ---
    PetscInt nLocal, nGlobal, nLocalMax;
    ierr = DMSwarmGetLocalSize(user->swarm, &nLocal); CHKERRQ(ierr);
    ierr = DMSwarmGetSize(user->swarm, &nGlobal); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&nLocal, &nLocalMax, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    PetscReal avg_per_rank = (size > 0) ? ((PetscReal)nGlobal / size) : 0.0;
    // Handle division by zero if there are no particles
    simCtx->particleLoadImbalance = (avg_per_rank > 1e-9) ? (nLocalMax / avg_per_rank) : 1.0;
 
 
    // --- 2. Number of Occupied Cells ---
    // This part requires access to the user->ParticleCount vector.
    PetscInt       local_occupied_cells = 0;
    PetscInt       global_occupied_cells;
    const PetscScalar *count_array;
    PetscInt       vec_local_size;
 
    ierr = VecGetLocalSize(user->ParticleCount, &vec_local_size); CHKERRQ(ierr);
    ierr = VecGetArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
 
    for (PetscInt i = 0; i < vec_local_size; ++i) {
        if (count_array[i] > 0.5) { // Use 0.5 to be safe with floating point
            local_occupied_cells++;
        }
    }
    ierr = VecRestoreArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
 
    ierr = MPI_Allreduce(&local_occupied_cells, &global_occupied_cells, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
    simCtx->occupiedCellCount = global_occupied_cells;
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Advanced Metrics: Imbalance=%.2f, OccupiedCells=%d\n", rank, simCtx->particleLoadImbalance, simCtx->occupiedCellCount);
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "LOG_PARTICLE_METRICS"
/**
 * @brief Implementation of \ref LOG_PARTICLE_METRICS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/logging.h`.
 * @see LOG_PARTICLE_METRICS()
 */
PetscErrorCode LOG_PARTICLE_METRICS(UserCtx *user, const char *stageName)
{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    SimCtx         *simCtx = user->simCtx;
    const char     *stage_label = (stageName && stageName[0] != '\0') ? stageName : "N/A";
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    PetscInt totalParticles;
    ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
 
    if (!rank) {
        FILE *f;
        char filen[PETSC_MAX_PATH_LEN + 64];
        ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Particle_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
        f = fopen(filen, "a");
        if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open particle log file: %s", filen);
 
        if (ftell(f) == 0) {
            PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10s | %-12s | %-10s | %-10s | %-10s | %-15s | %-10s | %-10s\n",
                            "Stage", "Timestep", "Total Ptls", "Lost", "Lost Total", "Migrated", "Occupied Cells", "Imbalance", "Mig Passes");
            PetscFPrintf(PETSC_COMM_SELF, f, "-------------------------------------------------------------------------------------------------------------------------------------------\n");
        }
        if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
            PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
        }
 
        PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10d | %-12d | %-10d | %-10d | %-10d | %-15d | %-10.2f | %-10d\n",
                        stage_label, (int)simCtx->step, (int)totalParticles, (int)simCtx->particlesLostLastStep,
                        (int)simCtx->particlesLostCumulative, (int)simCtx->particlesMigratedLastStep, (int)simCtx->occupiedCellCount,
                        (double)simCtx->particleLoadImbalance, (int)simCtx->migrationPassesLastStep);
        fclose(f);
    }
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PicurvOpenDiagnosticsCsv"
/**
 * @brief Implementation of \ref PicurvOpenDiagnosticsCsv().
 * @details Full API contract is documented with the header declaration in
 *          `include/logging.h`.
 */
PetscErrorCode PicurvOpenDiagnosticsCsv(const SimCtx *simCtx, const char *filename,
                                        const char *header, FILE **file)
{
    char  path[PETSC_MAX_PATH_LEN + 64];
    FILE *handle = NULL;
 
    PetscFunctionBeginUser;
    PetscCheck(simCtx && filename && header && file, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "A diagnostics CSV needs a context, a name, a header, and somewhere to "
               "return the handle.");
 
    PetscCall(PetscSNPrintf(path, sizeof(path), "%s/%s", simCtx->analysis_dir, filename));
    handle = fopen(path, "a");
    PetscCheck(handle != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
               "Unable to open diagnostics file '%s'.", path);
 
    if (ftell(handle) == 0) fprintf(handle, "%s\n", header);
    /* Written once, on the first step a continuation produces, so that a reader can tell
       a resumed run from a jump in the data itself. */
    if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
        fprintf(handle, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
    }
 
    *file = handle;
    PetscFunctionReturn(0);
}