variables.h

Go to the documentation of this file.

/**
 * @file variables.h
 * @author Vishal Kandala
 * @brief Main header file for a complex fluid dynamics solver.
 *
 * This file defines the core data structures, global variables, and function
 * prototypes for a CFD application. It combines features for a curvilinear,
 * multi-block, immersed boundary (IBM) solver with fluid-structure interaction (FSI),
 * a modern particle tracking system, a modular boundary condition framework,
 * and post-processing utilities. It heavily utilizes the PETSc library.
 *
 * *** REVISION NOTE ***
 * This version introduces a central `SimulationContext` struct to encapsulate all
 * simulation-wide configuration, replacing the legacy system of global variables.
 * The `UserCtx` struct has been refactored to hold only data specific to a
 * single grid level and now contains a back-pointer to the `SimulationContext`.
 * All struct definitions now follow the `typedef struct Name { ... } Name;` convention.
 */
 
#ifndef VARIABLES_H
#define VARIABLES_H
 
/* Build identity, generated by make so that a commit or dirty-state change rebuilds
 * the translation units that stamp it into binaries and checkpoints. The fallbacks
 * keep a build possible outside make; a binary reporting 0.0.0 is saying exactly that. */
#if defined(__has_include)
#  if __has_include("picurv_build_identity.h")
#    include "picurv_build_identity.h"
#  endif
#endif
#ifndef PICURV_RELEASE_VERSION
#define PICURV_RELEASE_VERSION "0.0.0"
#endif
#ifndef PICURV_GIT_COMMIT
#define PICURV_GIT_COMMIT "unknown"
#endif
#ifndef PICURV_BUILD_DIRTY
#define PICURV_BUILD_DIRTY "unknown"
#endif
 
/*================================================================================*
 *                        PETSC & SYSTEM LIBRARIES                                *
 *================================================================================*/
 
// --- PETSc Includes ---
#include "petscvec.h"
#include "petscdmda.h"
#include "petscksp.h"
#include "petscsnes.h"
#include "petscdmswarm.h"
 
// --- Standard C/C++ Includes ---
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
 
// --- C++ Specific Includes ---
#if defined(__cplusplus)
#include <vector>
#include <algorithm>
#include <cassert>
#endif
 
// Define a C-compatible interface for C++ compilers
#ifdef __cplusplus
extern "C" {
#endif
 
/*================================================================================*
 *                       MACROS & GLOBAL VARIABLES                                *
 *================================================================================*/
 
/// Coefficient controlling the temporal accuracy scheme (e.g., 1.5 for 2nd Order Backward Difference).
#define COEF_TIME_ACCURACY 1.5
 
/* --- NOTE: All legacy global 'extern' variables have been removed. --- */
/* They are now members of the SimulationContext struct. */
 
 
/*================================================================================*
 *                     CORE DATA STRUCTURES & TYPEDEFS                            *
 *================================================================================*/
 
//--------------------------------------------------------------------------------
//                 1. FORWARD DECLARATIONS & BASIC TYPES
//--------------------------------------------------------------------------------
 
// --- Forward Declarations ---
// These declarations allow pointers to these types before their full definition.
typedef struct SimCtx SimCtx;
struct ParticleFieldPlan;
struct PicurvWindow;
struct PicurvWindowStorage;
typedef struct UserCtx UserCtx;
typedef struct BC_Param_s BC_Param; /* Retains _s for linked list safety */
typedef struct BoundaryCondition BoundaryCondition;
typedef struct BoundaryFaceConfig BoundaryFaceConfig;
typedef struct IBMNodes IBMNodes;
typedef struct IBMVNodes IBMVNodes;
typedef struct FSInfo FSInfo;
typedef struct UserMG UserMG;
typedef struct BoundingBox BoundingBox;
typedef struct Cell Cell;
typedef struct Particle Particle;
typedef struct RankNeighbors RankNeighbors;
typedef struct RankCellInfo RankCellInfo;
typedef struct MigrationInfo MigrationInfo;
typedef struct Cstart Cstart;
typedef struct MGCtx MGCtx;
typedef struct PostProcessParams PostProcessParams;
typedef struct VTKFieldInfo VTKFieldInfo;
typedef struct VTKMetaData VTKMetaData;
typedef struct IBMInfo IBMInfo;
typedef struct SurfElmtInfo SurfElmtInfo;
typedef struct ScalingCtx ScalingCtx;
 
// --- Foundational Geometric and Data Types ---
 
/** @brief A 3D point or vector with PetscScalar components. */
typedef struct Cmpnts {
    PetscScalar x, y, z;
} Cmpnts;
 
/** @brief A 2D point or vector with PetscScalar components. */
typedef struct Cmpnts2 {
    PetscScalar x, y;
} Cmpnts2;
 
/** @brief A symmetric second-order tensor stored by its six independent components.
 *
 * The three missing components are recovered by symmetry: `yx == xy`, `zx == xz`,
 * and `zy == yz`. Contractions must therefore double the off-diagonal terms, which
 * is what `SymTensorContract()` exists to do in one place instead of at every call
 * site. Used for strain rates and for the Germano identity's Leonard and model
 * tensors.
 */
typedef struct SymTensor {
    PetscReal xx, xy, xz, yy, yz, zz;
} SymTensor;
 
/** @brief A 2D vector of PETSc real numbers (for geometry/coordinates). */
typedef struct Cpt2D {
  PetscReal x, y;
} Cpt2D;
 
/** @brief Boundary condition data used by the BC pipeline and ghost cell machinery. */
typedef struct BCS {
  Vec Ubcs;   ///< Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meaningful; interior is unused. Consumed by UpdateDummyCells for ghost cell extrapolation.
  Vec Uch;    ///< Characteristic velocity for boundary conditions.
} BCS;
 
 
//--------------------------------------------------------------------------------
//                 2. PARTICLE LOCATION SYSTEM ENUMS
//--------------------------------------------------------------------------------
 
/**
 * @enum  ParticleLocationStatus
 * @brief Defines the state of a particle with respect to its location and migration
 *        status during the iterative SettleParticles() process.
 */
typedef enum {
    NEEDS_LOCATION,
    ACTIVE_AND_LOCATED,
    MIGRATING_OUT,
    LOST,
    UNINITIALIZED
} ParticleLocationStatus;
 
/** @brief Enumerates the six faces of a cubic cell for distance calculations. */
typedef enum {
    LEFT = 0, RIGHT, BOTTOM, TOP, FRONT, BACK, NUM_FACES
} Face;
 
/** @brief Selects the algorithm used to populate a fresh Eulerian velocity field. */
typedef enum {
    IC_MODE_ZERO = 0,
    IC_MODE_CONSTANT_CARTESIAN = 1,
    IC_MODE_POISEUILLE = 2,
    IC_MODE_CONSTANT_STREAMWISE = 3,
    IC_MODE_FILE = 4
} InitialConditionMode;
 
/** @brief Selects the authoritative velocity representation in a staged file IC. */
typedef enum {
    IC_FIELD_UCAT = 0,
    IC_FIELD_UCONT = 1
} InitialConditionField;
 
 
//--------------------------------------------------------------------------------
//                 3. PARTICLE LOCATION SYSTEM STRUCTS
//--------------------------------------------------------------------------------
 
/** @brief Defines a 3D axis-aligned bounding box. */
typedef struct BoundingBox {
    Cmpnts min_coords; ///< Minimum x, y, z coordinates of the bounding box.
    Cmpnts max_coords; ///< Maximum x, y, z coordinates of the bounding box.
} BoundingBox;
 
/** @brief Defines the vertices of a single hexahedral grid cell. */
typedef struct Cell {
    Cmpnts vertices[8]; ///< Coordinates of the eight vertices of the cell.
} Cell;
 
/** @brief Defines a particle's core properties for Lagrangian tracking. */
typedef struct Particle {
    PetscInt64 PID;
    PetscInt cell[3];
    Cmpnts loc;
    Cmpnts vel;
    Cmpnts weights;
    ParticleLocationStatus location_status;
    PetscMPIInt destination_rank;
    PetscReal diffusivity;
    Cmpnts diffusivitygradient;
    PetscReal psi;
} Particle;
 
/** @brief Stores the MPI ranks of neighboring subdomains. */
typedef struct RankNeighbors {
    PetscMPIInt rank_xm, rank_xp;
    PetscMPIInt rank_ym, rank_yp;
    PetscMPIInt rank_zm, rank_zp;
} RankNeighbors;
 
/** @brief A lean struct to hold the global cell ownership range for a single MPI rank. */
typedef struct RankCellInfo {
    PetscInt xs_cell, ys_cell, zs_cell;
    PetscInt xm_cell, ym_cell, zm_cell;
} RankCellInfo;
 
/** @brief Information needed to migrate a single particle between MPI ranks. */
typedef struct MigrationInfo {
    PetscInt local_index;
    PetscInt target_rank;
} MigrationInfo;
 
/** @brief Verification-only analytical diffusivity override settings. */
typedef struct VerificationDiffusivityConfig {
    PetscBool enabled;
    char mode[32];
    char profile[32];
    PetscReal gamma0;
    PetscReal slope_x;
} VerificationDiffusivityConfig;
 
/** @brief Verification-only analytical scalar override settings. */
typedef struct VerificationScalarConfig {
    PetscBool enabled;
    char mode[32];
    char profile[32];
    PetscReal value;
    PetscReal phi0;
    PetscReal slope_x;
    PetscReal amplitude;
    PetscReal kx;
    PetscReal ky;
    PetscReal kz;
} VerificationScalarConfig;
 
/** @brief Aggregated per-timestep search instrumentation counters. */
typedef struct SearchMetricsState {
    PetscInt64 searchAttempts;
    PetscInt64 searchPopulation;
    PetscInt64 searchLocatedCount;
    PetscInt64 searchLostCount;
    PetscInt64 traversalStepsSum;
    PetscInt64 reSearchCount;
    PetscInt64 maxTraversalSteps;
    PetscInt64 maxTraversalFailCount;
    PetscInt64 tieBreakCount;
    PetscInt64 boundaryClampCount;
    PetscInt64 bboxGuessSuccessCount;
    PetscInt64 bboxGuessFallbackCount;
    PetscInt64 maxParticlePassDepth;
    PetscInt   currentSettlementPass;
} SearchMetricsState;
 
 
//--------------------------------------------------------------------------------
//                 4. BOUNDARY CONDITION SYSTEM ENUMS
//--------------------------------------------------------------------------------
 
/** @brief Identifies the six logical faces of a structured computational block. */
typedef enum {
    BC_FACE_NEG_X = 0, BC_FACE_POS_X = 1,
    BC_FACE_NEG_Y = 2, BC_FACE_POS_Y = 3,
    BC_FACE_NEG_Z = 4, BC_FACE_POS_Z = 5
} BCFace;
 
/** @brief Primary flow direction for streamwise IC and Poiseuille modes.
 *  Even values = positive-face-normal direction; odd = negative.
 *  axis = value / 2  (0=Xi, 1=Eta, 2=Zeta).
 *  sign = (value % 2 == 0) ? +1 : -1.
 */
typedef enum {
    FLOW_DIR_POS_XI   = 0,
    FLOW_DIR_NEG_XI   = 1,
    FLOW_DIR_POS_ETA  = 2,
    FLOW_DIR_NEG_ETA  = 3,
    FLOW_DIR_POS_ZETA = 4,
    FLOW_DIR_NEG_ZETA = 5,
    FLOW_DIR_UNSET    = -1
} FlowDirection;
 
/** @brief Defines the general mathematical/physical Category of a boundary. */
typedef enum {
    WALLFUNCTION       = -1,
    INTERFACE          = 0,
    WALL               = 1 , 
    MOVING_WALL        = 2,
    SYMMETRY           = 3,
    OUTLET             = 4,
    INLET              = 5, 
    FARFIELD           = 6,
    PERIODIC           = 7,
    CHARACTERISTIC_BC  = 8,
    ANALYTICAL_VORTEX  = 9,
    JUNCTION           = 10,
    ANNULUS            = 11,
    OGRID              = 12,
    RHEOLOGY           = 13,
    // Note:  Legacy 14 can be a JUNCTION with a specific handler.
} BCType;
 
/** @brief Defines the specific computational "strategy" for a boundary handler. */
typedef enum {
    BC_HANDLER_UNDEFINED = 0,
    BC_HANDLER_WALL_NOSLIP = 1, 
    BC_HANDLER_WALL_MOVING = 2,
    BC_HANDLER_SYMMETRY_PLANE = 3,
    BC_HANDLER_INLET_CONSTANT_VELOCITY = 4,
    BC_HANDLER_INLET_PARABOLIC = 5,
    BC_HANDLER_INLET_PROFILE_FROM_FILE = 6,
    BC_HANDLER_FARFIELD_NONREFLECTING = 9,
    BC_HANDLER_OUTLET_CONSERVATION = 10, 
    BC_HANDLER_OUTLET_PRESSURE = 11,
    BC_HANDLER_PERIODIC_GEOMETRIC = 12, 
    BC_HANDLER_INTERFACE_OVERSET = 13,
    BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX =  14,
    BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX = 15
  } BCHandlerType;
 
  
typedef enum{
    BC_PRIORITY_UNDEFINED = -1,
    BC_PRIORITY_INLET = 0,
    BC_PRIORITY_FARFIELD = 1,
    BC_PRIORITY_WALL = 2,
    BC_PRIORITY_OUTLET =3
} BCPriorityType;
 
//--------------------------------------------------------------------------------
//               5. BOUNDARY CONDITION SYSTEM STRUCTS
//--------------------------------------------------------------------------------
 
/** @brief A node in a linked list for storing key-value parameters from the bcs.dat file. */
struct BC_Param_s {
    char *key;
    char *value;
    struct BC_Param_s *next;
};
 
/** @brief Provides execution context for a boundary condition handler. */
typedef struct BCContext {
    UserCtx  *user;
    BCFace    face_id;
    const PetscReal *global_inflow_sum;
    const PetscReal *global_farfield_inflow_sum;
    const PetscReal *global_farfield_outflow_sum;
    const PetscReal *global_outflow_sum;
} BCContext;
 
/** @brief The "virtual table" struct for a boundary condition handler object. */
typedef struct BoundaryCondition {
    BCHandlerType type;
    BCPriorityType priority;
    void         *data;
    PetscErrorCode (*Initialize)(BoundaryCondition *self, BCContext *ctx);
    PetscErrorCode (*PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow);
    PetscErrorCode (*Apply)(BoundaryCondition *self, BCContext *ctx);
    PetscErrorCode (*PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow);
    PetscErrorCode (*UpdateUbcs)(BoundaryCondition *self,BCContext *ctx);
    PetscErrorCode (*Destroy)(BoundaryCondition *self);
} BoundaryCondition;
 
/** @brief Holds the complete configuration for one of the six boundary faces. */
typedef struct BoundaryFaceConfig {
    BCFace             face_id;
    BCType             mathematical_type;
    BCHandlerType      handler_type;
    BC_Param           *params;
    BoundaryCondition *handler;
} BoundaryFaceConfig;
 
//--------------------------------------------------------------------------------
//               6. IMBERSED BOUNDARY METHOD (IBM) & FSI STRUCTS
//--------------------------------------------------------------------------------
 
/** @brief Holds interpolation and distance information for a single IBM point. */
typedef struct IBMInfo {
  PetscInt  i1, j1, k1, i2, j2, k2, i3, j3, k3;
  PetscReal cr1, cr2, cr3;
  PetscReal d_i;
  PetscInt  imode;
  PetscInt  ni, nj, nk;
  PetscReal d_s;
  Cmpnts    pmin;
  PetscInt  cell;
  PetscReal cs1, cs2, cs3;
  PetscInt  i11, j11, k11, i22, j22, k22, i33, j33, k33;
  PetscReal cr11, cr22, cr33;
  PetscReal d_ii;
  PetscInt  iimode;
  PetscReal cs11, cs22, cs33;
  PetscInt  ii1, jj1, kk1, ii2, jj2, kk2, ii3, jj3, kk3;
  PetscReal ct1, ct2, ct3;
  PetscInt  smode;
  PetscInt  ii11, jj11, kk11, ii22, jj22, kk22, ii33, jj33, kk33;
  PetscReal ct11, ct22, ct33;
  PetscReal d_ss;
  PetscInt  ssmode;
} IBMInfo;
 
/** @brief Represents a collection of nodes forming a surface for the IBM. */
typedef struct IBMNodes {
  PetscInt  nbnumber;
  PetscInt  n_v, n_elmt;
  PetscInt  *nv1, *nv2, *nv3;
  PetscReal *nf_x, *nf_y, *nf_z;
  PetscReal *x_bp, *y_bp, *z_bp;
  PetscReal *x_bp0, *y_bp0, *z_bp0;
  PetscReal *x_bp_o, *y_bp_o, *z_bp_o;
  PetscReal *cf, *cfsum, *CCp;
  PetscReal x_bp_in[101][3270], y_bp_in[101][3270], z_bp_in[101][3270];
  Cmpnts    *u, *uold, *urm1;
  PetscReal *dA, *nt_x, *nt_y, *nt_z, *ns_x, *ns_y, *ns_z;
  PetscReal *cent_x, *cent_y, *cent_z;
  PetscReal *pres, *tau0, *tauN, *Bvel_u, *Bvel_v, *Bvel_w;
  PetscReal x_min, x_max, y_min, y_max, z_min, z_max;
  Cmpnts    *qvec;
  PetscReal *radvec;
} IBMNodes;
 
/** @brief Represents a collection of nodes forming a volume for the IBM. */
typedef struct IBMVNodes {
  PetscInt  nbnumber;
  PetscInt  n_v, n_elmt;
  PetscInt  *nv1, *nv2, *nv3, *nv4;
  PetscReal *x_bp, *y_bp, *z_bp;
  PetscReal *x_bp0, *y_bp0, *z_bp0;
  PetscReal *x_bp_o, *y_bp_o, *z_bp_o;
  Cmpnts    *u, *uold, *urm1;
  PetscReal V, *dV0;
  PetscReal *cent_x, *cent_y, *cent_z;
  PetscReal x_c, y_c, z_c;
  PetscReal J[3][3], I_inv[3][3];
} IBMVNodes;
 
/** @brief A generic C-style linked list node for integers. */
typedef struct node{
  PetscInt Node;
  struct node *next;
} node;
 
/** @brief Head of a generic C-style linked list. */
typedef struct list{
  node *head;
} List;
 
/** @brief A generic C-style linked list node for indices. */
typedef struct list_node {
  PetscInt  index;
  struct list_node *next;
} Node_List;
 
/** @brief A linked list node containing IBM interpolation info. */
typedef struct IBMListNode {
  IBMInfo ibm_intp;
  struct IBMListNode* next;
} IBMListNode;
 
/** @brief Head of a linked list for IBM data. */
typedef struct IBMList {
  IBMListNode *head;
} IBMList;
 
/** @brief Information about a surface element for FSI calculations. */
typedef struct SurfElmtInfo {
  PetscReal     P;
  PetscInt      n_P;
  PetscReal     Tow_ws, Tow_wt, Tow_wn;
  PetscInt      Clsnbpt_i, Clsnbpt_j, Clsnbpt_k;
  PetscInt      icell, jcell, kcell;
  PetscInt      FoundAroundcell, Need3rdPoint;
} SurfElmtInfo;
 
/** @brief Holds all data related to the state and motion of a body in FSI. */
typedef struct FSInfo {
  PetscReal    S_new[6], S_old[6], S_real[6], S_realm1[6];
  PetscReal    S_ang_n[6], S_ang_o[6], S_ang_r[6], S_ang_rm1[6];
  PetscReal    red_vel, damp, mu_s;
  PetscReal    F_x, F_y, F_z, A_tot;
  PetscReal    F_x_old, F_y_old, F_z_old;
  PetscReal    F_x_real, F_y_real, F_z_real;
  PetscReal    M_x, M_y, M_z;
  PetscReal    M_x_old, M_y_old, M_z_old;
  PetscReal    M_x_real, M_y_real, M_z_real;
  PetscReal    M_x_rm2, M_y_rm2, M_z_rm2;
  PetscReal    M_x_rm3, M_y_rm3, M_z_rm3;
  PetscReal    x_c, y_c, z_c;
  PetscReal    a_c[3];
  PetscReal    Mdpdn_x, Mdpdn_y, Mdpdn_z;
  PetscReal    Mdpdn_x_old, Mdpdn_y_old, Mdpdn_z_old;
  PetscReal    Power, clone;
  PetscInt     pbc[3];
  PetscReal    I_inv[3][3];
  PetscReal    L_n[3], L_o[3], L_r[3];
  PetscReal    alpha[3], acc[3];
  PetscReal    R[3][3], q[4], q_r[4];
  PetscReal    dS[6], dS_o[6], atk, atk_o;
  SurfElmtInfo *elmtinfo;
  IBMInfo      *fsi_intp;
  PetscReal    Max_xbc, Min_xbc, Max_ybc, Min_ybc, Max_zbc, Min_zbc;
  PetscInt     CV_ys, CV_ye, CV_zs, CV_ze;
} FSInfo;
 
/** @brief Defines prescribed body motion using splines. */
typedef struct Cstart {
  PetscInt   n_time, n_midp, n_subit;
  PetscReal  *x_midp, *y_midp, *z_midp;
  PetscReal  *x_com, *y_com, *head_ang;
  PetscReal  *s1, *s2, *s3, *st1, *st2, *st3;
  Mat        Mphi;
  PetscReal  xmin, xmax, ymin, ymax, zmin, zmax;
} Cstart;
 
//--------------------------------------------------------------------------------
//               7. LES TURBULENCE MODEL STRUCTS & ENUMS
//--------------------------------------------------------------------------------
/** @brief Identifies the subgrid-scale closure evaluated during a timestep.
 *
 * Add new enum values only when the parser, runtime dispatch, docs, and tests
 * are updated in the same change.
 */
typedef enum {
    NO_LES_MODEL = 0,
    CONSTANT_SMAGORINSKY = 1,
    DYNAMIC_SMAGORINSKY = 2,
    VREMAN = 3,
    WALE = 4
} LESModelType;
 
/** @brief Selects the wall model applied on WALL faces.
 *
 * The value doubles as the enable switch: ::WALL_FUNCTION_NONE disables the treatment,
 * and every other value both enables it and names the law. Configured independently of
 * the LES closure.
 *
 * Add new enum values only when the parser, runtime dispatch, docs, and tests are
 * updated in the same change.
 */
typedef enum {
    WALL_FUNCTION_NONE    = 0,
    WALL_FUNCTION_LOG_LAW = 1,
    WALL_FUNCTION_WERNER  = 2,
    WALL_FUNCTION_CABOT   = 3
} WallFunctionModel;
 
/** @brief Selects how a cell's grid filter width is derived from its metrics.
 *
 * `CUBE_ROOT_VOLUME` is exact for a cube and progressively underestimates the
 * width as a cell is stretched. The others use the cell's extents along its own grid
 * directions: `SCOTTI` corrects the cube root for the cell's aspect ratios, and
 * `MAX_EDGE` takes the longest extent. See @ref ComputeCellFilterWidth.
 */
typedef enum {
    LES_FILTER_WIDTH_CUBE_ROOT_VOLUME = 0,
    LES_FILTER_WIDTH_GEOMETRIC_MEAN   = 1,
    LES_FILTER_WIDTH_MAX_EDGE         = 2,
    LES_FILTER_WIDTH_SCOTTI           = 3
} LESFilterWidthModel;
 
/** @brief Selects the discrete test-filter kernel used by the dynamic procedure.
 *
 * `SIMPSON_IK` collapses onto the central eta-plane and is valid only when the
 * xi and zeta directions are homogeneous; the ingress layer enforces that.
 */
typedef enum {
    LES_TEST_FILTER_VOLUME_WEIGHTED_BOX = 0,
    LES_TEST_FILTER_SIMPSON_IK          = 1
} LESTestFilterKernel;
 
/** @brief Selects the set over which the Germano contractions are averaged.
 *
 * Lilly's least-squares closure assumes the model coefficient is constant across
 * the averaging set, so the set should span directions in which the flow really is
 * statistically homogeneous. `LOCAL` makes no such claim and is the only choice
 * valid on arbitrary geometry. `HOMOGENEOUS` derives its directions from the
 * periodic boundary pairs unless they are named explicitly.
 */
typedef enum {
    LES_AVERAGING_LOCAL       = 0,
    LES_AVERAGING_HOMOGENEOUS = 1,
    LES_AVERAGING_GLOBAL      = 2
} LESAveragingMode;
 
/** @brief Selects the admissible range imposed on the dynamic model coefficient.
 *
 * `NONE` keeps the coefficient signed, so subgrid-to-resolved energy transfer
 * (backscatter) survives; stability then rests on the total-viscosity floor
 * rather than on the coefficient's sign. The clipping modes discard the negative
 * tail, which removes backscatter and biases the mean dissipation upward.
 */
typedef enum {
    LES_CLIP_CLAMP         = 0,
    LES_CLIP_CLIP_NEGATIVE = 1,
    LES_CLIP_NONE          = 2
} LESClipMode;
 
/** @brief Every user-selectable parameter of the LES closure.
 *
 * Resolved once from the generated control file by @ref CreateSimulationContext and
 * read-only thereafter. The model selector itself stays on `SimCtx::les`, because
 * dispatch and field availability are decided from it throughout the runtime.
 */
typedef struct LESConfig {
    PetscInt            dynamic_frequency;       ///< Recompute the dynamic coefficient every N steps.
    PetscReal           constant_cs;             ///< Fixed Cs for CONSTANT_SMAGORINSKY; unused by the dynamic model.
    PetscReal           vreman_coefficient;      ///< Model constant c for VREMAN (Vreman 2004: 2.5 Cs^2).
    PetscReal           wale_coefficient;        ///< Model constant C_w for WALE (Nicoud & Ducros 1999).
    LESFilterWidthModel filter_width_model;      ///< How the grid filter width Delta is derived per cell.
    LESTestFilterKernel test_filter_kernel;      ///< Discrete test-filter stencil.
    PetscReal           test_filter_width_ratio; ///< Test-to-grid ratio per filtered direction; alpha is its square for the box, ratio^(4/3) for Simpson.
    LESAveragingMode    averaging_mode;          ///< Averaging set for the Germano contractions.
    PetscBool           averaging_direction[3];  ///< Averaged-over logical directions (xi, eta, zeta).
    LESClipMode         clip_mode;               ///< Admissible range for the coefficient.
    PetscReal           max_cs;                  ///< Ceiling on Cs under LES_CLIP_CLAMP.
    PetscReal           min_viscosity_ratio;     ///< Enforce nu + nu_t >= ratio * nu.
    PetscReal           yoshizawa_ci;            ///< Yoshizawa constant for the reported SGS kinetic energy.
    PetscBool           diagnostics_enabled;     ///< Append per-step coefficient statistics to the run log directory.
    PetscInt            diagnostics_cadence;     ///< Steps between diagnostic rows.
} LESConfig;
 
/** @brief Pre-clipping volume statistics captured by one dynamic-coefficient update.
 *
 * The backscattering and limited fractions describe the coefficient field as the
 * Germano contraction produced it. Clipping overwrites that state, and no stored
 * field preserves it, so it is recorded here on the way past and consumed by
 * @ref LogLESDiagnostics.
 */
typedef struct LESDiagnosticsState {
    PetscReal contraction_lm;      ///< Volume-weighted sum of `L_ij M_ij` over fluid cells.
    PetscReal contraction_mm;      ///< Volume-weighted sum of `M_ij M_ij` over fluid cells.
    PetscReal backscatter_volume;  ///< Fluid volume whose raw coefficient was negative.
    PetscReal limited_volume;      ///< Fluid volume whose coefficient the clip modified.
    PetscReal fluid_volume;        ///< Total fluid volume sampled.
    PetscBool valid;               ///< Set once a dynamic update has populated this state.
} LESDiagnosticsState;
 
/** @brief Near-wall statistics captured by one wall-model pass.
 *
 * The wall model evaluates the wall distance and the friction velocity together at each
 * corrected cell. `y+` is the number that says whether the model is being asked to work
 * in the range it is valid for, and it cannot be recovered afterwards from the corrected
 * velocity field, which no longer carries the law that produced it. Recomputing the wall
 * distance in the diagnostic would duplicate geometry the correction already has in hand,
 * so it is accumulated on the way past and consumed by @ref LogWallModelDiagnostics.
 */
typedef struct WallModelDiagnosticsState {
    PetscReal friction_velocity_sum;  ///< Sum of `u_tau` over corrected cells.
    PetscReal friction_velocity_sq;   ///< Sum of `u_tau^2`, for the RMS.
    PetscReal friction_velocity_min;  ///< Smallest `u_tau` this rank corrected.
    PetscReal friction_velocity_max;  ///< Largest `u_tau` this rank corrected.
    PetscReal wall_distance_sum;      ///< Sum of the first-cell wall distance.
    PetscReal y_plus_sum;             ///< Sum of `u_tau * y / nu` over corrected cells.
    PetscReal y_plus_max;             ///< Largest first-cell `y+` this rank corrected.
    PetscReal wall_viscosity_sum;     ///< Sum of the effective wall eddy viscosity.
    PetscInt  cells;                  ///< Number of cells this rank corrected.
} WallModelDiagnosticsState;
//--------------------------------------------------------------------------------
//               8. MULTIGRID, SOLVERS AND POST-PROCESSING STRUCTS AND ENUMS
//--------------------------------------------------------------------------------
 
/** @brief Enumerator to identify the implemented momentum solver strategies.
 *
 * Add new enum values only when the parser, runtime dispatch, docs, and tests
 * are updated in the same change.
 */
typedef enum{
    MOMENTUM_SOLVER_EXPLICIT_RK = 0,
    MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK = 1,
    MOMENTUM_SOLVER_NEWTON_KRYLOV = 2,
    /** @deprecated Use MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK. */
    MOMENTUM_SOLVER_DUALTIME_PICARD_RK4 = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK
} MomentumSolverType;
 
/** @brief Selects the runtime solution-convergence diagnostics mode. */
typedef enum{
    SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC = 0,
    SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC = 1,
    SOLUTION_CONVERGENCE_STATISTICAL_STEADY = 2,
    SOLUTION_CONVERGENCE_TRANSIENT = 3
} SolutionConvergenceMode;
 
/** @brief Enumerator to identify the particle initialization strategy. */
typedef enum {
    PARTICLE_INIT_SURFACE_RANDOM = 0, /**< Random placement on the inlet face */
    PARTICLE_INIT_VOLUME         = 1, /**< Random volumetric distribution across the domain */
    PARTICLE_INIT_POINT_SOURCE   = 2, /**< All particles at a fixed (psrc_x,psrc_y,psrc_z) — for validation */
    PARTICLE_INIT_SURFACE_EDGES  = 3  /**< Deterministic placement at inlet face edges */
} ParticleInitializationType;
 
/**
 * @brief Selects the grid-to-particle interpolation method.
 *
 * - INTERP_TRILINEAR: Direct trilinear from 8 nearest cell centers (second-order on curvilinear grids).
 * - INTERP_CORNER_AVERAGED: Legacy center-to-corner average, then trilinear from corners.
 */
typedef enum {
    INTERP_TRILINEAR       = 0,
    INTERP_CORNER_AVERAGED = 1
} InterpolationMethod;
 
/** @brief Context for Multigrid operations. */
typedef struct MGCtx {
  UserCtx  *user;
  PetscInt thislevel;
  DM       packer;
} MGCtx;
 
/** @brief User-level context for managing the entire multigrid hierarchy. */
typedef struct UserMG {
  PetscInt  mglevels;
  PetscInt  thislevel;
  PetscBool isc, jsc, ksc;
  MGCtx     *mgctx;
  DM        packer;
  SNES      snespacker;
} UserMG;
 
#define MAX_PIPELINE_LENGTH 1024
#define MAX_FIELD_LIST_LENGTH 1024
#define MAX_FILENAME_LENGTH 256
/// Defines the maximum number of data fields for VTK point data.
#define MAX_POINT_DATA_FIELDS 20 
#define MAX_VTK_FIELD_NAME_LENGTH 64 ///< Maximum length for VTK field names.
/**
 * @brief Holds all configuration parameters for a post-processing run.
 *        This is an enhanced version combining command-line and file-based settings.
 */
typedef struct PostProcessParams {
    // -- Source Directory --- (For Data)
    char source_dir[PETSC_MAX_PATH_LEN];
  
    // --- Time Controls (can be set by command line or file) ---
    PetscInt startTime;
    PetscInt endTime;
    PetscInt timeStep;
    /** File listing the exact steps to process, one per line; empty uses the time controls. */
    char step_list_file[PETSC_MAX_PATH_LEN];
    PetscBool outputParticles;
 
    // --- Configuration primarily from the .cfg file ---
    char process_pipeline[MAX_PIPELINE_LENGTH];
    char output_fields_instantaneous[MAX_FIELD_LIST_LENGTH];
    char output_prefix[MAX_FILENAME_LENGTH];
    char particle_pipeline[MAX_PIPELINE_LENGTH];
    char particle_fields[MAX_FIELD_LIST_LENGTH];
    char particle_output_prefix[MAX_FILENAME_LENGTH];
    PetscInt particle_output_freq;
 
    // --- Statistics Pipeline (global reductions → CSV files) ---
    char statistics_pipeline[MAX_PIPELINE_LENGTH];      /**< e.g. "ComputeMSD" */
    char statistics_output_prefix[MAX_FILENAME_LENGTH]; /**< basename for CSV output, e.g. "Stats" */
 
    // --- Field Statistics (accumulated window state → derived fields) ---
    /** Comma-separated window names to derive; empty disables the pipeline. */
    char field_statistics_windows[MAX_FIELD_LIST_LENGTH];
    /** Comma-separated outputs: mean, reynolds_stress, rms, tke, flux. */
    char field_statistics_outputs[MAX_FIELD_LIST_LENGTH];
    /** Comma-separated formats: vtk for derived fields, csv for the convergence history. */
    char field_statistics_formats[MAX_FIELD_LIST_LENGTH];
    /** Analysis prefix for CSV summaries; VTK continues to use output_prefix. */
    char field_statistics_output_prefix[PETSC_MAX_PATH_LEN];
    /** Whether derived output leaves non-dimensional form, from global_operations.dimensionalize. */
    PetscBool dimensionalize;
    /** Committed step supplying the state; negative means the step being processed. */
    PetscInt field_statistics_source_step;
 
    // --- Legacy settings ---
    char eulerianExt[8]; // from original PostProcessParams (repurposed for PreCheckAndResize() as the input file extension.)
    char particleExt[8]; // from original PostProcessParams
 
    // --- Processing Parameters ---
    PetscInt reference[3]; // Reference point for normalizing any field against.
    
} PostProcessParams;
 
/**
 * @brief Stores all necessary information for a single data array in a VTK file.
 */
typedef struct VTKFieldInfo {
    char        name[MAX_VTK_FIELD_NAME_LENGTH]; // Name of the data field
    PetscInt    num_components; // 1 for scalar, 3 for vector
    PetscScalar* data;          // Pointer to the raw data array (must be freed by caller)
} VTKFieldInfo;
 
/** @brief Enumerates the type of VTK file to be written. */
typedef enum {
    VTK_STRUCTURED,
    VTK_POLYDATA
} VTKFileType;
 
typedef struct VTKMetaData {
    VTKFileType     fileType;
    PetscInt        mx, my, mz;
    PetscInt        npoints;
    PetscScalar*    coords;
    VTKFieldInfo    point_data_fields[MAX_POINT_DATA_FIELDS];
    PetscInt        num_point_data_fields;
    PetscInt*       connectivity;
    PetscInt*       offsets;
} VTKMetaData;
 
/**
 * @brief Defines the execution mode of the application.
 */
typedef enum {
    EXEC_MODE_SOLVER,          // The application is running as the main fluid solver.
    EXEC_MODE_POSTPROCESSOR,   // The application is running as the post-processor.
    EXEC_MODE_UNKNOWN          // Default/error state.
} ExecutionMode;
 
//-------------------------------------------------------------------------------
//             9. SCALING AND DIMENSIONAL ANALYSIS STRUCTS
//-------------------------------------------------------------------------------
typedef struct ScalingCtx{
  PetscReal L_ref;
  PetscReal U_ref;
  PetscReal rho_ref;
  PetscReal P_ref;
}ScalingCtx;
//-------------------------------------------------------------------------------
/*================================================================================*
 *                        MAIN APPLICATION CONTEXT                                *
 *================================================================================*/
 
/**
 * @brief The master context for the entire simulation.
 *
 * This struct encapsulates all global configuration flags, physical parameters,
 * simulation control settings, and top-level data objects. It replaces the
 * numerous global variables from the legacy codebase, providing a single,
 * explicit source of truth for the simulation's setup.
 */
typedef struct SimCtx {
 
    //================ Group 1: Parallelism & MPI Information ================
    PetscMPIInt rank;
    PetscMPIInt size;
 
    //================ Group 2: Simulation Control, Time, and I/O ================
 
    PetscInt  step;
    PetscReal ti;
    PetscInt  StartStep;
    PetscInt  StepsToRun;
    PetscInt  tiout;
    PetscInt  particleConsoleOutputFreq;
    PetscReal StartTime;
    PetscReal dt;
    PetscBool OnlySetup;
    PetscBool continueMode;
    PetscViewer logviewer;
    ExecutionMode exec_mode;
    char eulerianSource[PETSC_MAX_PATH_LEN];
    char restart_dir[PETSC_MAX_PATH_LEN];
    char output_dir[PETSC_MAX_PATH_LEN];
    char log_dir[PETSC_MAX_PATH_LEN];
    char analysis_dir[PETSC_MAX_PATH_LEN];
    char _io_context_buffer[PETSC_MAX_PATH_LEN]; // Persistent store for I/O context strings.
    char *current_io_directory; // Pointer into the above buffer.
    char checkpointGeometrySHA256[65];
    PetscBool checkpointGeometryHashReady;
    PetscBool restartHistoryAvailable;
 
    //================ Group 3: High-Level Physics & Model Selection Flags ================
    PetscInt  immersed, movefsi, rotatefsi;   ///< Refused at setup: immersed boundaries and moving bodies are not implemented.
    PetscInt  invicid, TwoD, moveframe, rotateframe; ///< moveframe/rotateframe are refused at setup.
    char AnalyticalSolutionType[PETSC_MAX_PATH_LEN];
  
    //================ Group 4: Immersed-Body Flux Corrections ================
    PetscInt  MHV, LV;                         ///< Heart-valve flux corrections for immersed bodies; refused at setup.
 
    //================ Group 5: Solver & Numerics Parameters ================
    MomentumSolverType mom_solver_type;
    PetscInt  mom_max_pseudo_steps;
    PetscReal mom_atol, mom_rtol, mom_resid_atol, mom_resid_rtol, imp_stol;
    PetscInt  mglevels, mg_preItr, mg_poItr;
    PetscInt  poisson;
    PetscInt  central;
 
    PetscReal ren, pseudo_cfl, cdisx, cdisy, cdisz;
    PetscReal pseudo_cfl_reduction_factor, pseudo_cfl_growth_factor; // st, vnn
    PetscReal max_pseudo_cfl, min_pseudo_cfl; // New addition for adaptive pseudo-CFL
    PetscReal mom_dt_jameson_residual_norm_noise_allowance_factor; // New addition for divergence detection
    PetscBool no_pseudo_cfl_backtrack; /* Diagnostic: commit all finite RK trials, skip ratio-based rollback */
    PetscReal mom_ratio_ema_alpha;    /* EMA smoothing coefficient for trial-ratio rejection (0=no smooth, 1=raw) */
    PetscBool mom_last_converged;     /* Set after each momentum solve: true if pseudo-time converged */
    PetscReal mom_last_lambda_max;   /* Global spectral radius [1/s] from the last momentum solve; for diagnostics */
    PetscBool mom_nk_monitor_history; /* Write PICurv's structured per-Newton-iteration residual history. */
    PetscBool ps_ksp_pic_monitor_true_residual; // Parsed once from options for custom Poisson monitor logging.
    InitialConditionMode initialConditionMode;
    InitialConditionField initialConditionField;
    char initialConditionDirectory[PETSC_MAX_PATH_LEN];
    Cmpnts    InitialConstantContra;  /* Cartesian (u,v,w) for constant-cartesian IC */
    FlowDirection flowDirection;      /* explicit flow direction for streamwise/Poiseuille */
    PetscReal icVelocityPhysical;     /* scalar speed for curvilinear Constant and Poiseuille */
    Cmpnts    AnalyticalUniformVelocity;
    PetscBool solutionConvergenceEnabled;
    SolutionConvergenceMode solutionConvergenceMode;
    PetscInt  solutionConvergencePeriodSteps;
    PetscInt  solutionConvergenceWindowSteps;
    PetscInt  solutionConvergenceSamplesRecorded;
    PetscReal *solutionConvergenceMeanSpeedHistory;
    PetscReal *solutionConvergenceMeanKEHistory;
    /* --- Field Statistics (window lifecycle; accumulator binding arrives with ingress) --- */
    PetscBool      fieldStatisticsEnabled;
    PetscInt       fieldStatisticsWindowCount;
    struct PicurvWindow *fieldStatisticsWindows;
    PetscInt       statisticsConsoleOutputFreq;
    /* Whether saved window state is resumed rather than restarted from zero.
     * Defaults to the run's continue mode; --restart-from requires an explicit
     * request because it may branch onto a different physical trajectory. */
    PetscBool      fieldStatisticsContinue;
    VerificationDiffusivityConfig verificationDiffusivity;
    VerificationScalarConfig verificationScalar;
    
    //================ Group 6: Physical & Geometric Parameters ================
    PetscReal Flux_in, angle,max_angle;
    PetscReal CMx_c, CMy_c, CMz_c;
    PetscReal psrc_x, psrc_y, psrc_z;  /**< Point source location for PARTICLE_INIT_POINT_SOURCE */
    ScalingCtx scaling;
    PetscReal  wall_roughness_height;
    PetscReal schmidt_number, Turbulent_schmidt_number;
    PetscReal iem_constant;            /**< IEM mixing constant C_IEM in Omega = C_IEM Gamma / Delta^2 (default 2.0). */
 
    //================ Group 7: Grid, Domain, and Boundary Condition Settings ================
    PetscInt  block_number, inletprofile;
    PetscInt  i_periodic, j_periodic, k_periodic, pseudo_periodic;
    PetscBool generate_grid;        
    char grid_file[PETSC_MAX_PATH_LEN];
    PetscInt da_procs_x, da_procs_y, da_procs_z;
    PetscInt num_bcs_files;
    char **bcs_files;
    PetscReal  FluxInSum, FluxOutSum,Fluxsum,FarFluxInSum,FarFluxOutSum;
    // Turbulent Flat Channel with Forcing term
    PetscReal drivingForceMagnitude,forceScalingFactor;
    /* Physical step whose smoothed driving force `drivingForceMagnitude` holds
     * (-1 = none yet). ComputeDrivenChannelFlowSource() runs once per ComputeRHS,
     * i.e. per RK stage or per Newton residual evaluation, so its EMA must be
     * advanced against this rather than on every call. */
    PetscInt  drivingForceStep;
    PetscReal targetVolumetricFlux; // for DRIVEN flow.
    /* True once a driven-flow target is fixed for the run. `constant_flux` sets
     * it at initialization from the bcs file; `initial_flux` sets it at the
     * first PreStep from the starting field. Checkpointed so a restart holds the
     * original target instead of re-measuring a drifted one. */
    PetscBool drivenFluxTargetLatched;
    PetscReal bulkVelocityCorrection;
    PetscReal boundaryVelocityCorrection;
    /* The planar-average volumetric flux and cross-section area the controller set
     * `bulkVelocityCorrection` from this step. Kept so the driven-flow diagnostics can
     * report them without repeating the controller's global reductions. */
    PetscReal drivenFluxMeasured, drivenFluxArea;
    PetscReal  AreaInSum, AreaOutSum;
    PetscInt   ccc;
    PetscReal  ratio;
  
    //================ Group 8: Turbulence Modeling (LES) ================
    PetscInt  les;                  ///< Active LES closure; an ::LESModelType value.
    PetscInt  wallfunction;         ///< Enable wall functions on WALL faces.
    PetscInt  les_gradient_model;   ///< Add the Clark gradient (tensor-diffusivity) term to the viscous flux.
    LESConfig les_config;           ///< Parameters of the LES closure selected by `les`.
  
    //================ Group 9: Particle / DMSwarm Data & Settings ================
    PetscInt  np;
    PetscInt  particleRandomSeed; ///< Base seed for every particle RNG stream (-particle_random_seed).
    DM        dm_swarm;
    BoundingBox *bboxlist;
    ParticleInitializationType ParticleInitialization;
    InterpolationMethod interpolationMethod;
    char particleRestartMode[16];
    PetscInt particlesLostLastStep;
    PetscInt particlesLostCumulative;
    PetscReal particlesLostScalarLastStep; /**< Sum of Psi over the particles removed this step. */
    struct ParticleFieldPlan *particleFieldPlan; /**< Configured initial particle-field values, or NULL. */
    PetscInt migrationPassesLastStep;
    PetscInt particlesMigratedLastStep;
    PetscInt occupiedCellCount;
    PetscReal particleLoadImbalance;
    SearchMetricsState searchMetrics;
    PetscRandom BrownianMotionRNG;
    PetscReal C_IEM;
 
    //================ Group 10: Immersed Boundary & FSI Data Object Pointers ================
    IBMNodes  *ibm;
    IBMVNodes *ibmv;
    FSInfo    *fsi;
 
    //================ Group 11: Top-Level Managers, Custom Configuration,Logging/Monitoring ====
    UserMG    usermg;
    char      allowedFile[PETSC_MAX_PATH_LEN];
    PetscBool useCfg;
    char      **allowedFuncs;
    PetscInt  nAllowed;
    PetscInt  LoggingFrequency;
    /* Net volume flux into the domain carried by the momentum step's velocity before the
     * pressure correction (inflow minus outflow, through every boundary face): the domain
     * integral of the pressure equation's source. With Neumann pressure boundaries it must
     * be near zero for that equation to have a solution. Unlike the continuity log's
     * Net Flux, which sums only the inlet and outlet faces, it covers all boundaries. */
    PetscReal poissonSourceImbalance;
    PetscReal MaxDiv;
    PetscInt  MaxDivFlatArg, MaxDivx,MaxDivy,MaxDivz;
    // Profiling 
    char      profilingSelectedFuncsFile[PETSC_MAX_PATH_LEN];
    PetscBool useProfilingSelectedFuncsCfg;
    char      **profilingSelectedFuncs;
    PetscInt  nProfilingSelectedFuncs;
    char      profilingTimestepMode[32];
    char      profilingTimestepFile[PETSC_MAX_PATH_LEN];
    PetscBool profilingFinalSummary;
    PetscBool walltimeGuardEnabled;
    PetscBool walltimeGuardActive;
    PetscInt  walltimeGuardWarmupSteps;
    PetscReal walltimeGuardMultiplier;
    PetscReal walltimeGuardMinSeconds;
    PetscReal walltimeGuardEstimatorAlpha;
    PetscReal walltimeGuardJobStartEpochSeconds;
    PetscReal walltimeGuardLimitSeconds;
    PetscInt  walltimeGuardCompletedSteps;
    PetscReal walltimeGuardWarmupTotalSeconds;
    PetscReal walltimeGuardWarmupAverageSeconds;
    PetscBool walltimeGuardHasEWMA;
    PetscReal walltimeGuardEWMASeconds;
    PetscReal walltimeGuardLatestStepSeconds;
    PetscBool runtimeMemoryLogEnabled;               /**< Enable the rank-reduced runtime memory log. */
    char      runtimeMemoryLogFile[PETSC_MAX_PATH_LEN]; /**< File name written under log_dir. */
    PetscBool runtimeMemoryLogStarted;               /**< True after rank 0 writes the log header. */
    PetscBool runtimeMemoryLogHasPrevious;           /**< True after the first process-memory sample. */
    PetscReal runtimeMemoryLogPreviousProcessMB;     /**< Previous local process memory sample in MB. */
    //================ Group 12: Post-Processing =================================================
    char      PostprocessingControlFile[PETSC_MAX_PATH_LEN];
    PostProcessParams *pps;
    
   //=============== Group 13: Miscellaneous =============================================
   PetscReal    r[101], tin[101], uinr[101][1001];
 
} SimCtx;
 
/** @deprecated Use mom_dt_jameson_residual_norm_noise_allowance_factor. */
#define mom_dt_rk4_residual_norm_noise_allowance_factor mom_dt_jameson_residual_norm_noise_allowance_factor
 
 
/**
 * @brief User-defined context containing data specific to a single computational
 *        grid level. This is the primary data structure passed to low-level
 *        numerical and solver functions.
 */
typedef struct UserCtx {
 
    // --- The Critical Link to Global Configuration ---
    SimCtx *simCtx;  ///< Back-pointer to the master simulation context.
 
    // --- Grid, Geometry & Parallelization (Per-Level) ---
    DM da, fda;
    /* Symmetric second-order tensor DM (dof 6), mirroring da's decomposition.
     * Field-statistics products need six components at each point; splitting them
     * across six scalar vectors would cost six memory streams in the per-step
     * accumulation loop and six collective gathers per checkpoint. */
    DM fda6;
    DMDALocalInfo info;
    AO ao;
    PetscInt IM, JM, KM; // No.of grid points.
    PetscReal Max_X, Max_Y, Max_Z, Min_X, Min_Y, Min_Z;
    BoundingBox bbox;
    RankNeighbors neighbors;
    PetscInt GridOrientation, isc, jsc, ksc, _this;
    PetscReal rx, ry, rz;
    PetscInt cgrid;
    Cmpnts periodic_translation[3];
    PetscBool periodic_translation_valid[3];
 
    // --- Boundary Conditions & Interfacing ---
    BoundaryFaceConfig boundary_faces[6];
    PetscBool inletFaceDefined;
    BCFace    identifiedInletBCFace;
    BCS       Bcs;
    /* Wall-model friction velocity. Global is authoritative - the wall model writes
       it alongside the global Ucat it corrects - and the local view is its ghosted
       image, refreshed through the field catalog like every other cell field. */
    Vec       Friction_Velocity, lFriction_Velocity;
    /* Eddy viscosity the wall model needs at its wall face for the discrete viscous
       flux to deliver the stress the model computed. Derived, rebuilt by every wall
       pass, and never restored from a checkpoint. */
    Vec       Nu_Wall, lNu_Wall;
 
    // --- Primary Flow Fields (Global & Local Views) ---
    Vec Ucont, lUcont, Ucat, lUcat, P, lP, Phi, lPhi, Nvert, lNvert;
 
    // --- Secondary Flow & Derived Fields ---
    Vec Diffusivity, lDiffusivity;
    Vec DiffusivityGradient, lDiffusivityGradient;
 
    // --- Time-Stepping & Solver Workspace Fields ---
    Vec Ucont_o, lUcont_o, Ucat_o, P_o, Nvert_o, lNvert_o;
    Vec Ucont_rm1, lUcont_rm1, Rhs, dUcont, pUcont;
    Vec *solutionConvergencePeriodicUcatRef;
    Vec *solutionConvergencePeriodicPRef;
    /* Corner-staging workspace, one pair per block size, so neither the DM nor
     * the degree of freedom has to be rediscovered from a cached vector. */
    Vec CellScalarAtCorner, lCellScalarAtCorner;   /* da-based,  dof 1 */
    Vec CellVectorAtCorner, lCellVectorAtCorner;   /* fda-based, dof 3 */
    /* Post-processing staging vectors. A derived statistic is config-counted and has
     * no compile-time offset of its own, so it is written here to be addressed by
     * name through the existing ghost and nodal-average paths. The nodal companions
     * are reused field by field, because the VTK writer copies each one out before
     * the next is produced. */
    Vec PostScalar, lPostScalar, PostScalarNodal;  /* da-based,  dof 1 */
    Vec PostVector, lPostVector, PostVectorNodal;  /* fda-based, dof 3 */
    /* Per-window statistics accumulators, one entry per configured window. */
    struct PicurvWindowStorage *fieldStatisticsStorage;
 
    // --- Pressure-Poisson System ---
    /* Built once by PoissonSolver_Multigrid() and kept for the run. A, nullsp and R live
       on every level (R below the finest); MR and MP are the transfers between this level
       and the next coarser one; ksp and B are held on the finest level only. */
    Mat A; KSP ksp; MatNullSpace nullsp;
    Vec B,R;
    Mat MR, MP;
  
  
    // --- Grid Metrics (Global and Local) ---
  Vec Cent, lCent, Csi, Eta, Zet, Aj, lCsi, lEta, lZet, lAj, GridSpace,lGridSpace;
  Vec   Centx,Centy,Centz;
  Vec   lCentx,lCenty,lCentz;
  Vec   ICsi, IEta, IZet, IAj, lICsi, lIEta, lIZet, lIAj;
  Vec   JCsi, JEta, JZet, JAj, lJCsi, lJEta, lJZet, lJAj;
  Vec   KCsi, KEta, KZet, KAj, lKCsi, lKEta, lKZet, lKAj;
 
  // --- Turbulence Modeling (LES) ---
  Vec Nu_t, lNu_t, CS, lCs, Distance;
  LESDiagnosticsState les_diagnostics; ///< Pre-clipping statistics from the last dynamic update.
  WallModelDiagnosticsState wall_diagnostics; ///< Near-wall statistics from the last wall-model pass.
  PetscInt wall_yplus_excursions; ///< Consecutive diagnostic samples with the first cell outside the selected law's valid y+ range.
 
  // --- Immersed Boundary Method (IBM) ---
  IBMNodes *ibm; IBMList *ibmlist;
 
  // --- Multigrid Hierarchy ---
  PetscInt thislevel, mglevels;
  UserCtx *user_f, *user_c;
 
  // --- Particle System ---
  DM swarm;
  RankCellInfo *RankCellInfoMap;
  Vec ParticleCount, lParticleCount;
  Vec Psi, lPsi; //scalar dummy to demonstrate scatter.
  
  // --- Post-Processing ---
  DM  post_swarm;
  Vec P_nodal;
  Vec Ucat_nodal;
  Vec Qcrit, lQcrit;   ///< Cell-centred Q-criterion and its ghosted copy for nodal averaging.
  Vec Qcrit_nodal;     ///< Q-criterion averaged to grid nodes; the field a .vts can place correctly.
  Vec Psi_nodal;
  
} UserCtx;
 
#ifdef __cplusplus
}
#endif
 
#endif // VARIABLES_H