grid.c

Go to the documentation of this file.

// in src/grid.c
 
#include "grid.h"
#include "statistics_window.h"
#include "logging.h"
 
#define BBOX_TOLERANCE 1e-6
 
#undef __FUNCT__
#define __FUNCT__ "ParseAndSetGridInputs"
/**
 * @brief Parse grid-generation options and store the resulting geometry settings in the context.
 */
static PetscErrorCode ParseAndSetGridInputs(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx; // Get the global context via the back-pointer
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
    if(strcmp(simCtx->eulerianSource,"analytical")==0 &&
       AnalyticalTypeRequiresCustomGeometry(simCtx->AnalyticalSolutionType)){
        ierr = SetAnalyticalGridInfo(user); CHKERRQ(ierr);
    } else if (simCtx->generate_grid) {
        if (strcmp(simCtx->eulerianSource, "analytical") == 0) {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                           "Rank %d: Analytical type '%s' uses programmatic grid ingestion for block %d.\n",
                           simCtx->rank, simCtx->AnalyticalSolutionType, user->_this);
        } else {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: Block %d is programmatically generated. Calling generation parser.\n", simCtx->rank, user->_this);
        }
        ierr = ReadGridGenerationInputs(user); CHKERRQ(ierr);
    } else {
        if (strcmp(simCtx->eulerianSource, "analytical") == 0) {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                           "Rank %d: Analytical type '%s' uses file-based grid ingestion for block %d.\n",
                           simCtx->rank, simCtx->AnalyticalSolutionType, user->_this);
        } else {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: Block %d is file-based. Calling file parser.\n", simCtx->rank, user->_this);
        }
        ierr = ReadGridFile(user); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "DefineAllGridDimensions"
/**
 * @brief Internal helper implementation: `DefineAllGridDimensions()`.
 * @details Local to this translation unit.
 */
PetscErrorCode DefineAllGridDimensions(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    PetscInt       nblk = simCtx->block_number;
    UserCtx        *finest_users;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (simCtx->usermg.mglevels == 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "MG levels not set. Cannot get finest_users.");
    }
    // Get the UserCtx array for the finest grid level
    finest_users = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Defining grid dimensions for %d blocks...\n", nblk);
    if (strcmp(simCtx->eulerianSource, "analytical") == 0 &&
        AnalyticalTypeRequiresCustomGeometry(simCtx->AnalyticalSolutionType)) {
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "Analytical type '%s' requires custom geometry; preloading finest-grid IM/JM/KM once.\n",
                  simCtx->AnalyticalSolutionType);
        ierr = PopulateFinestUserGridResolutionFromOptions(finest_users, nblk); CHKERRQ(ierr);
    }
 
    // Loop over each block to configure its grid dimensions and geometry.
    for (PetscInt bi = 0; bi < nblk; bi++) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: --- Configuring Geometry for Block %d ---\n", simCtx->rank, bi);
 
        // Before calling any helpers, set the block index in the context.
        // This makes the UserCtx self-aware of which block it represents.
        LOG_ALLOW(GLOBAL,LOG_DEBUG,"finest_users->_this = %d, bi = %d\n",finest_users[bi]._this,bi);
        //finest_user[bi]._this = bi;
 
        // Call the helper function for this specific block. It can now derive
        // all necessary information from the UserCtx pointer it receives.
        ierr = ParseAndSetGridInputs(&finest_users[bi]); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "CreateCompatibleBlockDM"
/**
 * @brief Implementation of \ref CreateCompatibleBlockDM().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see CreateCompatibleBlockDM()
 */
PetscErrorCode CreateCompatibleBlockDM(DM source, PetscInt dof, DM *result)
{
    PetscErrorCode ierr;
    PetscInt       dim = 0, M = 0, N = 0, P = 0, m = 0, n = 0, p = 0, source_dof = 0, stencil_width = 0;
    DMBoundaryType bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;
    DMDAStencilType stencil_type = DMDA_STENCIL_BOX;
    const PetscInt *lx = NULL, *ly = NULL, *lz = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(source != NULL && result != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Source DM and output pointer are required.");
    PetscCheck(dof > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "A compatible DM needs a positive degree of freedom, got %" PetscInt_FMT ".", dof);
 
    ierr = DMDAGetInfo(source, &dim, &M, &N, &P, &m, &n, &p, &source_dof, &stencil_width,
                       &bx, &by, &bz, &stencil_type); CHKERRQ(ierr);
    ierr = DMDAGetOwnershipRanges(source, &lx, &ly, &lz); CHKERRQ(ierr);
    ierr = DMDACreate3d(PetscObjectComm((PetscObject)source), bx, by, bz, stencil_type,
                        M, N, P, m, n, p, dof, stencil_width, lx, ly, lz, result); CHKERRQ(ierr);
    ierr = DMSetUp(*result); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG,
                   "Created a dof-%d DM mirroring the block decomposition (%dx%dx%d ranks).\n",
                   (int)dof, (int)m, (int)n, (int)p);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InitializeSingleGridDM"
/**
 * @brief Create and configure one PETSc DMDA for a multigrid level.
 */
static PetscErrorCode InitializeSingleGridDM(UserCtx *user, UserCtx *coarse_user)
{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
 
    DMBoundaryType xperiod = (simCtx->i_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    DMBoundaryType yperiod = (simCtx->j_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    DMBoundaryType zperiod = (simCtx->k_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    PetscInt stencil_width = (simCtx->i_periodic || simCtx->j_periodic || simCtx->k_periodic) ? 3:2; // Stencil width is 2 in the legacy code
 
    PetscInt *lx = NULL, *ly = NULL, *lz = NULL;
    PetscInt m, n, p;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (coarse_user) {
        // --- This is a FINE grid; it must be aligned with the COARSE grid ---
        LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: [Aligning DM] for block %d level %d (size %dx%dx%d) with level %d\n", simCtx->rank, user->_this, user->thislevel, user->IM, user->JM, user->KM, coarse_user->thislevel);
 
        DMDAGetInfo(coarse_user->da, NULL, NULL, NULL, NULL, &m, &n, &p, NULL, NULL, NULL, NULL, NULL, NULL);
        LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d:   Coarse grid processor decomposition is %d x %d x %d\n", simCtx->rank, m, n, p);
 
        // This is the core logic from MGDACreate to ensure processor alignment.
        PetscInt *lx_contrib, *ly_contrib, *lz_contrib;
        ierr = PetscMalloc3(m, &lx_contrib, n, &ly_contrib, p, &lz_contrib); CHKERRQ(ierr);
        ierr = PetscMemzero(lx_contrib, m * sizeof(PetscInt)); CHKERRQ(ierr);
        ierr = PetscMemzero(ly_contrib, n * sizeof(PetscInt)); CHKERRQ(ierr);
        ierr = PetscMemzero(lz_contrib, p * sizeof(PetscInt)); CHKERRQ(ierr);
 
        DMDALocalInfo info;
        DMDAGetLocalInfo(coarse_user->da, &info);
        PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
        PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
        PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
        
        PetscMPIInt rank;
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
        PetscInt proc_i = rank % m;
        PetscInt proc_j = (rank / m) % n;
        PetscInt proc_k = rank / (m * n);
 
        // --- X-Direction Logic (Identical to MGDACreate) ---
        if (user->isc) lx_contrib[proc_i] = (xe - xs);
        else {
            if (m == 1) lx_contrib[0] = user->IM + 1;
            else if (xs == 0) lx_contrib[0] = 2 * xe - 1;
            else if (xe == mx) lx_contrib[proc_i] = user->IM + 1 - (2 * xs - 1);
            else lx_contrib[proc_i] = (xe - xs) * 2;
        }
 
        // --- Y-Direction Logic (Identical to MGDACreate) ---
        if (user->jsc) ly_contrib[proc_j] = (ye - ys);
        else {
            if (n == 1) ly_contrib[0] = user->JM + 1;
            else if (ys == 0) ly_contrib[0] = 2 * ye - 1;
            else if (ye == my) ly_contrib[proc_j] = user->JM + 1 - (2 * ys - 1);
            else ly_contrib[proc_j] = (ye - ys) * 2;
        }
 
        // --- Z-Direction Logic (Identical to MGDACreate) ---
        if (user->ksc) lz_contrib[proc_k] = (ze - zs);
        else {
            if (p == 1) lz_contrib[0] = user->KM + 1;
            else if (zs == 0) lz_contrib[0] = 2 * ze - 1;
            else if (ze == mz) lz_contrib[proc_k] = user->KM + 1 - (2 * zs - 1);
            else lz_contrib[proc_k] = (ze - zs) * 2;
        }
        LOG_ALLOW_SYNC(LOCAL, LOG_VERBOSE, "Rank %d:   Calculated this rank's node contribution to fine grid: lx=%d, ly=%d, lz=%d\n", simCtx->rank, lx_contrib[proc_i], ly_contrib[proc_j], lz_contrib[proc_k]);
 
        // Allocate the final distribution arrays and Allreduce to get the global distribution
        ierr = PetscMalloc3(m, &lx, n, &ly, p, &lz); CHKERRQ(ierr);
        ierr = MPI_Allreduce(lx_contrib, lx, m, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(ly_contrib, ly, n, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(lz_contrib, lz, p, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        ierr = PetscFree3(lx_contrib, ly_contrib, lz_contrib); CHKERRQ(ierr);
 
    } else {
        // --- CASE 2: This is the COARSEST grid; use default or user-specified decomposition ---
        if(simCtx->exec_mode == EXEC_MODE_SOLVER){
 
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating coarsest DM for block %d level %d (size %dx%dx%d)\n", simCtx->rank, user->_this, user->thislevel, user->IM, user->JM, user->KM);
            m = simCtx->da_procs_x;
            n = simCtx->da_procs_y;
            p = simCtx->da_procs_z;
        
        } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
 
          LOG_ALLOW(GLOBAL,LOG_ERROR,"Currently Only Single Rank is supported. \n");
          
          m = n = p = PETSC_DECIDE;
            
        } 
        // lx, ly, lz are NULL, so DMDACreate3d will use the m,n,p values.
    }
 
    // --- Create the DMDA for the current UserCtx ---
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d:   Calling DMDACreate3d...\n", simCtx->rank);
    ierr = DMDACreate3d(PETSC_COMM_WORLD, xperiod, yperiod, zperiod, DMDA_STENCIL_BOX,
                          user->IM + 1, user->JM + 1, user->KM + 1,
                          m, n, p,
                          1, stencil_width, lx, ly, lz, &user->da); CHKERRQ(ierr);
    
    if (coarse_user) {
        ierr = PetscFree3(lx, ly, lz); CHKERRQ(ierr);
    }
    
    // --- Standard DM setup applicable to all levels ---
    ierr = DMSetUp(user->da); CHKERRQ(ierr);
    ierr = DMGetCoordinateDM(user->da, &user->fda); CHKERRQ(ierr);
    /* The symmetric-tensor DM exists only to carry statistics products, and only the
     * finest level accumulates them. Configuration is resolved before grid setup, so
     * a run without statistics never builds it. */
    if (FieldStatisticsIsActive(simCtx) && user->thislevel == simCtx->mglevels - 1) {
        ierr = CreateCompatibleBlockDM(user->da, 6, &user->fda6); CHKERRQ(ierr);
    }
    ierr = DMDASetUniformCoordinates(user->da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(user->da, &user->info); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d:   DM creation for block %d level %d complete.\n", simCtx->rank, user->_this, user->thislevel);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "InitializeAllGridDMs"
/**
 * @brief Internal helper implementation: `InitializeAllGridDMs()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InitializeAllGridDMs(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx = usermg->mgctx;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "Pre-scanning BCs to identify domain periodicity.\n");
    ierr = DeterminePeriodicity(simCtx); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Creating DMDA objects for all levels and blocks...\n");
 
    // --- Part 1: Calculate Coarse Grid Dimensions & VALIDATE ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating and validating coarse grid dimensions...\n");
    for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *user_coarse = &mgctx[level].user[bi];
            UserCtx *user_fine   = &mgctx[level + 1].user[bi];
 
            user_coarse->IM = user_fine->isc ? user_fine->IM : (user_fine->IM + 1) / 2;
            user_coarse->JM = user_fine->jsc ? user_fine->JM : (user_fine->JM + 1) / 2;
            user_coarse->KM = user_fine->ksc ? user_fine->KM : (user_fine->KM + 1) / 2;
 
            LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Block %d, Level %d dims calculated: %d x %d x %d\n",
                      simCtx->rank, bi, level, user_coarse->IM, user_coarse->JM, user_coarse->KM);
 
            // Validation check from legacy MGDACreate to ensure coarsening is possible
            PetscInt check_i = user_coarse->IM * (2 - user_coarse->isc) - (user_fine->IM + 1 - user_coarse->isc);
            PetscInt check_j = user_coarse->JM * (2 - user_coarse->jsc) - (user_fine->JM + 1 - user_coarse->jsc);
            PetscInt check_k = user_coarse->KM * (2 - user_coarse->ksc) - (user_fine->KM + 1 - user_coarse->ksc);
 
            if (check_i + check_j + check_k != 0) {
          //      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
          //           "Grid at level %d, block %d cannot be coarsened from %dx%dx%d to %dx%dx%d with the given semi-coarsening flags. Check grid dimensions.",
              //          level, bi, user_fine->IM, user_fine->JM, user_fine->KM, user_coarse->IM, user_coarse->JM, user_coarse->KM);
          LOG(GLOBAL,LOG_WARNING,"WARNING: Grid at level %d, block %d can't be consistently coarsened further.\n", level, bi);
            }
        }
    }
 
    // --- Part 2: Create DMs from Coarse to Fine for each Block ---
    for (PetscInt bi = 0; bi < nblk; bi++) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "--- Creating DMs for Block %d ---\n", bi);
 
        // Create the coarsest level DM first (passing NULL for the coarse_user)
        ierr = InitializeSingleGridDM(&mgctx[0].user[bi], NULL); CHKERRQ(ierr);
 
        // Create finer level DMs, passing the next-coarser context for alignment
        for (PetscInt level = 1; level < usermg->mglevels; level++) {
            ierr = InitializeSingleGridDM(&mgctx[level].user[bi], &mgctx[level-1].user[bi]); CHKERRQ(ierr);
        }
    }
 
    // --- Optional: View the finest DM for debugging verification ---
    if (get_log_level() >= LOG_DEBUG) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "--- Viewing Finest DMDA (Level %d, Block 0) ---\n", usermg->mglevels - 1);
        ierr = DMView(mgctx[usermg->mglevels - 1].user[0].da, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "DMDA object creation complete.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
// Forward declarations for the static helper functions within this file.
static PetscErrorCode SetFinestLevelCoordinates(UserCtx *user);
static PetscErrorCode GenerateAndSetCoordinates(UserCtx *user);
static PetscErrorCode ReadAndSetCoordinates(UserCtx *user, FILE *fd);
static PetscErrorCode RestrictCoordinates(UserCtx *coarse_user, UserCtx *fine_user);
 
#undef __FUNCT__
#define __FUNCT__ "AssignAllGridCoordinates"
/**
 * @brief Internal helper implementation: `AssignAllGridCoordinates()`.
 * @details Local to this translation unit.
 */
PetscErrorCode AssignAllGridCoordinates(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Assigning physical coordinates to all grid DMs...\n");
 
    // --- Part 1: Populate the Finest Grid Level ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Setting coordinates for the finest grid level (%d)...\n", usermg->mglevels - 1);
    for (PetscInt bi = 0; bi < nblk; bi++) {
        UserCtx *fine_user = &usermg->mgctx[usermg->mglevels - 1].user[bi];
        ierr = SetFinestLevelCoordinates(fine_user); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL,LOG_TRACE,"The Finest level coordinates for block %d have been set.\n",bi);
        if(get_log_level()==LOG_VERBOSE && is_function_allowed(__FUNCT__)==true){
            ierr = LOG_FIELD_MIN_MAX(fine_user, FIELD_ID_COORDINATES);
        }
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finest level coordinates have been set for all blocks.\n");
 
    // --- Part 2: Restrict Coordinates to Coarser Levels ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Restricting coordinates to coarser grid levels...\n");
    for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *coarse_user = &usermg->mgctx[level].user[bi];
            UserCtx *fine_user   = &usermg->mgctx[level + 1].user[bi];
            ierr = RestrictCoordinates(coarse_user, fine_user); CHKERRQ(ierr);
 
            LOG_ALLOW(GLOBAL,LOG_TRACE,"Coordinates restricted to block %d level %d.\n",bi,level);
            if(get_log_level()==LOG_VERBOSE && is_function_allowed(__FUNCT__)==true) {
                ierr = LOG_FIELD_MIN_MAX(coarse_user, FIELD_ID_COORDINATES);
            }
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Physical coordinates assigned to all grid levels and blocks.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Returns one Cartesian component from a coordinate/vector value.
 */
static inline PetscReal CoordinateComponent(Cmpnts value, PetscInt component)
{
    if (component == 0) return value.x;
    if (component == 1) return value.y;
    return value.z;
}
 
#undef __FUNCT__
#define __FUNCT__ "ValidatePeriodicGeometry"
/**
 * @brief Implementation of \ref ValidatePeriodicGeometry().
 * @details Full API contract is documented with the header declaration in
 *          `include/grid.h`.
 */
PetscErrorCode ValidatePeriodicGeometry(UserCtx *user)
{
    const BCFace neg_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace pos_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
    const char axis_names[3] = {'X', 'Y', 'Z'};
    const Cmpnts ***coor = NULL;
    Vec lcoor = NULL;
    DMDALocalInfo info;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(DMDAGetLocalInfo(user->da, &info));
    PetscCall(DMGetCoordinatesLocal(user->da, &lcoor));
    PetscCheck(lcoor != NULL, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "Cannot validate periodic geometry before local coordinates are assigned.");
    for (PetscInt axis = 0; axis < 3; axis++) {
        const PetscBool neg_periodic =
            user->boundary_faces[neg_faces[axis]].mathematical_type == PERIODIC;
        const PetscBool pos_periodic =
            user->boundary_faces[pos_faces[axis]].mathematical_type == PERIODIC;
        PetscReal local_min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
        PetscReal local_max[3] = {-PETSC_MAX_REAL, -PETSC_MAX_REAL, -PETSC_MAX_REAL};
        PetscReal global_min[3], global_max[3];
        PetscInt local_count = 0, global_count = 0;
 
        user->periodic_translation_valid[axis] = PETSC_FALSE;
        user->periodic_translation[axis] = (Cmpnts){0.0, 0.0, 0.0};
        if (!neg_periodic && !pos_periodic) continue;
 
        PetscCheck(neg_periodic && pos_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "Periodic geometry in the %c direction requires paired negative and positive faces.",
                   axis_names[axis]);
        const PetscInt axis_size = axis == 0 ? info.mx : (axis == 1 ? info.my : info.mz);
        PetscCheck(axis_size >= 5, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "%c-periodic geometry on block %d level %d requires at least four physical "
                   "nodes in that direction; found %d.",
                   axis_names[axis], user->_this, user->thislevel, axis_size - 1);
 
        PetscCall(DMDAVecGetArrayRead(user->fda, lcoor, &coor));
        if (axis == 0 && info.xs == 0) {
            for (PetscInt k = PetscMax(info.zs, 0); k < PetscMin(info.zs + info.zm, info.mz - 1); k++) {
                for (PetscInt j = PetscMax(info.ys, 0); j < PetscMin(info.ys + info.ym, info.my - 1); j++) {
                    const Cmpnts delta = {
                        coor[k][j][-2].x - coor[k][j][0].x,
                        coor[k][j][-2].y - coor[k][j][0].y,
                        coor[k][j][-2].z - coor[k][j][0].z
                    };
                    for (PetscInt c = 0; c < 3; c++) {
                        local_min[c] = PetscMin(local_min[c], CoordinateComponent(delta, c));
                        local_max[c] = PetscMax(local_max[c], CoordinateComponent(delta, c));
                    }
                    local_count++;
                }
            }
        } else if (axis == 1 && info.ys == 0) {
            for (PetscInt k = PetscMax(info.zs, 0); k < PetscMin(info.zs + info.zm, info.mz - 1); k++) {
                for (PetscInt i = PetscMax(info.xs, 0); i < PetscMin(info.xs + info.xm, info.mx - 1); i++) {
                    const Cmpnts delta = {
                        coor[k][-2][i].x - coor[k][0][i].x,
                        coor[k][-2][i].y - coor[k][0][i].y,
                        coor[k][-2][i].z - coor[k][0][i].z
                    };
                    for (PetscInt c = 0; c < 3; c++) {
                        local_min[c] = PetscMin(local_min[c], CoordinateComponent(delta, c));
                        local_max[c] = PetscMax(local_max[c], CoordinateComponent(delta, c));
                    }
                    local_count++;
                }
            }
        } else if (axis == 2 && info.zs == 0) {
            for (PetscInt j = PetscMax(info.ys, 0); j < PetscMin(info.ys + info.ym, info.my - 1); j++) {
                for (PetscInt i = PetscMax(info.xs, 0); i < PetscMin(info.xs + info.xm, info.mx - 1); i++) {
                    const Cmpnts delta = {
                        coor[-2][j][i].x - coor[0][j][i].x,
                        coor[-2][j][i].y - coor[0][j][i].y,
                        coor[-2][j][i].z - coor[0][j][i].z
                    };
                    for (PetscInt c = 0; c < 3; c++) {
                        local_min[c] = PetscMin(local_min[c], CoordinateComponent(delta, c));
                        local_max[c] = PetscMax(local_max[c], CoordinateComponent(delta, c));
                    }
                    local_count++;
                }
            }
        }
        PetscCall(DMDAVecRestoreArrayRead(user->fda, lcoor, &coor));
 
        PetscCallMPI(MPI_Allreduce(local_min, global_min, 3, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD));
        PetscCallMPI(MPI_Allreduce(local_max, global_max, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
        PetscCallMPI(MPI_Allreduce(&local_count, &global_count, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
        PetscCheck(global_count > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                   "No physical seam nodes were available to validate %c-periodic geometry.",
                   axis_names[axis]);
 
        PetscReal translation[3];
        PetscReal scale = 1.0;
        PetscReal max_mismatch = 0.0;
        for (PetscInt c = 0; c < 3; c++) {
            translation[c] = 0.5 * (global_min[c] + global_max[c]);
            scale = PetscMax(scale, PetscAbsReal(translation[c]));
            max_mismatch = PetscMax(max_mismatch, global_max[c] - global_min[c]);
        }
        const PetscReal tolerance = 1.0e-9 * scale + 100.0 * PETSC_MACHINE_EPSILON;
        const PetscReal magnitude = PetscSqrtReal(
            PetscSqr(translation[0]) + PetscSqr(translation[1]) + PetscSqr(translation[2]));
 
        PetscCheck(max_mismatch <= tolerance, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "Unsupported %c-periodic geometry on block %d level %d: opposite physical "
                   "surfaces are not related by one constant translation. Maximum component "
                   "mismatch is %.12e (tolerance %.12e).",
                   axis_names[axis], user->_this, user->thislevel,
                   (double)max_mismatch, (double)tolerance);
        PetscCheck(magnitude > tolerance, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "Unsupported %c-periodic geometry on block %d level %d: seam translation "
                   "magnitude %.12e is zero or too small.",
                   axis_names[axis], user->_this, user->thislevel, (double)magnitude);
 
        user->periodic_translation[axis] =
            (Cmpnts){translation[0], translation[1], translation[2]};
        user->periodic_translation_valid[axis] = PETSC_TRUE;
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "Validated %c-periodic geometry for block %d level %d with translation "
                  "(%.12e, %.12e, %.12e).\n",
                  axis_names[axis], user->_this, user->thislevel,
                  (double)translation[0], (double)translation[1], (double)translation[2]);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "SetFinestLevelCoordinates"
/**
 * @brief Attach the generated physical coordinates to the finest-level DMDA.
 */
static PetscErrorCode SetFinestLevelCoordinates(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting finest level coordinates for block %d...\n", simCtx->rank, user->_this);
 
    if (simCtx->generate_grid) {
        ierr = GenerateAndSetCoordinates(user); CHKERRQ(ierr);
    } else {
 
      FILE *grid_file_handle = NULL;
      // Only Rank 0 opens the file.
      if (simCtx->rank == 0) {
    grid_file_handle = fopen(simCtx->grid_file, "r");
    if (!grid_file_handle) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file: %s", simCtx->grid_file);
        
    // Now, on Rank 0, we skip the entire header section once.
    // This is the logic from your modern code's AssignGridCoordinates.
    PetscInt headerLines = simCtx->block_number + 2; // 1 for nblk, plus one for each block's dims
    char dummy_buffer[2048];
    for (PetscInt s = 0; s < headerLines; ++s) {
      if (!fgets(dummy_buffer, sizeof(dummy_buffer), grid_file_handle)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unexpected EOF while skipping grid header");
      }
    }
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank 0: Skipped %d header lines, now at coordinate data.\n", headerLines);
      }
 
      // We now call the coordinate reader, passing the file handle.
      // It's responsible for reading its block's data and broadcasting.
      ierr = ReadAndSetCoordinates(user, grid_file_handle); CHKERRQ(ierr);
        
      // Only Rank 0, which opened the file, should close it.
      if (simCtx->rank == 0) {
    fclose(grid_file_handle);
      }
    }
    
    // Populate local ghost coordinates from the owned global coordinates.
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Scattering coordinates to update ghost nodes for block %d...\n", simCtx->rank, user->_this);
    ierr = UpdateLocalGhosts(user, FIELD_ID_COORDINATES); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
/**
 * @brief Map a uniform logical coordinate to its configured stretched physical coordinate.
 */
static inline PetscReal ComputeStretchedCoord(PetscInt i, PetscInt N, PetscReal L, PetscReal r)
{
    if (N <=1) return 0.0;
    PetscReal fraction = (PetscReal)i / ((PetscReal)N - 1.0);
    if (PetscAbsReal(r - 1.0) < 1.0e-9) { // Use a tolerance for float comparison
        return L * fraction;
    } else {
        return L * (PetscPowReal(r, fraction) - 1.0) / (r - 1.0);
    }
}
 
#undef __FUNCT__
#define __FUNCT__ "GenerateAndSetCoordinates"
/**
 * @brief Generate analytic grid coordinates and install them on the finest grid level.
 */
static PetscErrorCode GenerateAndSetCoordinates(UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Cmpnts       ***coor;
    Vec            gCoor;
    
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Generating coordinates for block %d...\n", user->simCtx->rank, user->_this);
    
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMGetCoordinates(user->da, &gCoor); CHKERRQ(ierr);
 
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Local Info for block %d - X range - [%d,%d], Y range - [%d,%d], Z range - [%d,%d]\n",
              user->simCtx->rank, user->_this, xs, xe, ys, ye, zs, ze);
    
    LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Local Info for block %d - X domain - [%.4f,%.4f], Y range - [%.4f,%.4f], Z range - [%.4f,%.4f]\n",
              user->simCtx->rank, user->_this, user->Min_X,user->Max_X,user->Min_Y,user->Max_Y,user->Min_Z,user->Max_Z);          
              
    ierr = VecSet(gCoor, 0.0); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, gCoor, &coor); CHKERRQ(ierr);
    
    PetscReal Lx = user->Max_X - user->Min_X;
    PetscReal Ly = user->Max_Y - user->Min_Y;
    PetscReal Lz = user->Max_Z - user->Min_Z;
 
    // Loop over the nodes owned by this process.
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                if(k < user->KM && j < user->JM && i < user->IM){
                    coor[k][j][i].x = user->Min_X + ComputeStretchedCoord(i, user->IM, Lx, user->rx);
                    coor[k][j][i].y = user->Min_Y + ComputeStretchedCoord(j, user->JM, Ly, user->ry);
                    coor[k][j][i].z = user->Min_Z + ComputeStretchedCoord(k, user->KM, Lz, user->rz);
                }
            }  
        }
    }
 
    /// DEBUG: This verifies the presence of a last "unphysical" layer of coordinates.
    /*
    PetscInt KM = user->KM;
    for (PetscInt j = ys; j < ye; j++){
        for(PetscInt i = xs; i < xe; i++){
            LOG_ALLOW(GLOBAL,LOG_DEBUG,"coor[%d][%d][%d].(x,y,z) = %le,%le,%le",KM,j,i,coor[KM][j][i].x,coor[KM][j][i].y,coor[KM][j][i].z);
        }
    }
    */
 
 
 
    ierr = DMDAVecRestoreArray(user->fda, gCoor, &coor); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ReadAndSetCoordinates"
/**
 * @brief Load grid coordinates from the input file and install them on the finest level.
 */
static PetscErrorCode ReadAndSetCoordinates(UserCtx *user, FILE *fd)
{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
    PetscMPIInt    rank = simCtx->rank;
    PetscInt       block_index = user->_this;
    PetscInt       IM = user->IM, JM = user->JM, KM = user->KM;
    DMDALocalInfo  info;
    Cmpnts       ***coor;
    Vec            gCoor;
    PetscReal      *gc = NULL; // Global coordinate buffer, allocated on all ranks
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Reading interleaved coordinates from file for block %d...\n",
              simCtx->rank, block_index);
 
    // 1. Allocate the buffer on ALL ranks to receive the broadcast data.
    // PetscInt n_nodes = (IM + 1) * (JM + 1) * (KM + 1);
    PetscInt n_nodes = (IM) * (JM) * (KM);
    ierr = PetscMalloc1(3 * n_nodes, &gc); CHKERRQ(ierr);
 
    // 2. Only Rank 0 opens the file and reads the data.
    if (rank == 0) {
        if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Recieved a NULL file handle.\n");
                
        // Read the coordinate data for the CURRENT block.
        for (PetscInt k = 0; k < KM; k++) {
            for (PetscInt j = 0; j < JM; j++) {
                for (PetscInt i = 0; i < IM; i++) {
                    PetscInt base_index = 3 * ((k * (JM) + j) * (IM) + i);
                    if (fscanf(fd, "%le %le %le\n", &gc[base_index], &gc[base_index + 1], &gc[base_index + 2]) != 3) {
                        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Error reading coordinates for node (i,j,k)=(%d,%d,%d) in block %d", i, j, k, block_index);
                    }
                }
            }
        }
        
    }
 
    // 3. Broadcast the coordinate block for the current block to all other processes.
    ierr = MPI_Bcast(gc, 3 * n_nodes, MPIU_REAL, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // 4. Each rank populates its owned portion of the global coordinate vector.
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMGetCoordinates(user->da, &gCoor); CHKERRQ(ierr);
    ierr = VecSet(gCoor, 0.0); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, gCoor, &coor); CHKERRQ(ierr);
 
    for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
            for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
          if(k< KM && j < JM && i < IM){
                PetscInt base_idx = 3 * ((k * (JM) + j) * (IM) + i);
                coor[k][j][i].x = gc[base_idx];
                coor[k][j][i].y = gc[base_idx + 1];
                coor[k][j][i].z = gc[base_idx + 2];
          }
        }
        }
    }
    
    // 5. Clean up and restore.
    ierr = DMDAVecRestoreArray(user->fda, gCoor, &coor); CHKERRQ(ierr);
    ierr = PetscFree(gc); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "RestrictCoordinates"
/**
 * @brief Restrict finest-grid coordinates onto each coarser multigrid level.
 */
static PetscErrorCode RestrictCoordinates(UserCtx *coarse_user, UserCtx *fine_user)
{
    PetscErrorCode ierr;
    Vec            c_gCoor, f_lCoor;
    Cmpnts       ***c_coor;
    const Cmpnts ***f_coor; // Use const for read-only access
    DMDALocalInfo  c_info;
    PetscInt       ih, jh, kh; // Fine-grid indices corresponding to coarse-grid i,j,k
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Restricting coords from level %d to level %d for block %d\n",
                   fine_user->simCtx->rank, fine_user->thislevel, coarse_user->thislevel, coarse_user->_this);
 
    ierr = DMDAGetLocalInfo(coarse_user->da, &c_info); CHKERRQ(ierr);
 
    ierr = DMGetCoordinates(coarse_user->da, &c_gCoor); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(fine_user->da, &f_lCoor); CHKERRQ(ierr);
    
    ierr = VecSet(c_gCoor, 0.0); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(coarse_user->fda, c_gCoor, &c_coor); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(fine_user->fda, f_lCoor, &f_coor); CHKERRQ(ierr);
 
    // Get the local owned range of the coarse grid.
    PetscInt xs = c_info.xs, xe = c_info.xs + c_info.xm;
    PetscInt ys = c_info.ys, ye = c_info.ys + c_info.ym;
    PetscInt zs = c_info.zs, ze = c_info.zs + c_info.zm;
 
    // Get the global dimensions of the coarse grid.
    PetscInt mx = c_info.mx, my = c_info.my, mz = c_info.mz;
 
    // If this process owns the maximum boundary node, contract the loop by one
    // to prevent the index doubling `2*i` from going out of bounds.
    // This is also ensuring we do not manipulate the unphysical layer of coors present in the finest level.
    if (xe == mx) xe--;
    if (ye == my) ye--;
    if (ze == mz) ze--;
    
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                // Determine the corresponding parent node index on the FINE grid,
                // respecting the semi-coarsening flags of the FINE grid's UserCtx.
                ih = coarse_user->isc ? i : 2 * i;
                jh = coarse_user->jsc ? j : 2 * j;
        kh = coarse_user->ksc ? k : 2 * k;
 
        //  LOG_ALLOW(GLOBAL,LOG_DEBUG," [kh][ih][jh] = %d,%d,%d  - k,j,i = %d,%d,%d.\n",kh,jh,ih,k,j,i);
 
                c_coor[k][j][i] = f_coor[kh][jh][ih];
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(coarse_user->fda, c_gCoor, &c_coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(fine_user->fda, f_lCoor, &f_coor); CHKERRQ(ierr);
 
    // Populate the coarse-grid local ghost coordinates from the restricted global coordinates.
    ierr = UpdateLocalGhosts(coarse_user, FIELD_ID_COORDINATES); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ComputeLocalBoundingBox"
/**
 * @brief Implementation of \ref ComputeLocalBoundingBox().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see ComputeLocalBoundingBox()
 */
PetscErrorCode ComputeLocalBoundingBox(UserCtx *user, BoundingBox *localBBox)
{
    PetscErrorCode ierr;
    PetscInt i, j, k;
    PetscMPIInt rank;
    PetscInt xs, ys, zs, xe, ye, ze;
    DMDALocalInfo info;
    Vec coordinates;
    Cmpnts ***coordArray;
    Cmpnts minCoords, maxCoords;
 
    PetscFunctionBeginUser;
    
    PROFILE_FUNCTION_BEGIN;
 
    // Start of function execution
    LOG_ALLOW(GLOBAL, LOG_INFO, "Entering the function.\n");
 
    // Validate input Pointers
    if (!user) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Input 'user' Pointer is NULL.\n");
        PROFILE_FUNCTION_END;
        return PETSC_ERR_ARG_NULL;
    }
    if (!localBBox) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Output 'localBBox' Pointer is NULL.\n");
        PROFILE_FUNCTION_END;
        return PETSC_ERR_ARG_NULL;
    }
 
    // Get MPI rank
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // Get the local coordinates vector from the DMDA
    ierr = DMGetCoordinatesLocal(user->da, &coordinates);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error getting local coordinates vector.\n");
        PROFILE_FUNCTION_END;
        return ierr;
    }
 
    if (!coordinates) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Coordinates vector is NULL.\n");
        PROFILE_FUNCTION_END;
        return PETSC_ERR_ARG_NULL;
    }
 
    // Access the coordinate array for reading
    ierr = DMDAVecGetArrayRead(user->fda, coordinates, &coordArray);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error accessing coordinate array.\n");
        PROFILE_FUNCTION_END;
        return ierr;
    }
 
    // Get the local grid information (indices and sizes)
    ierr = DMDAGetLocalInfo(user->da, &info);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error getting DMDA local info.\n");
        PROFILE_FUNCTION_END;
        return ierr;
    }
 
    
    xs = info.gxs; xe = xs + info.gxm;
    ys = info.gys; ye = ys + info.gym;
    zs = info.gzs; ze = zs + info.gzm;
    
    /*
    xs = info.xs; xe = xs + info.xm;
    ys = info.ys; ye = ys + info.ym;
    zs = info.zs; ze = zs + info.zm;
    */
    
    // Initialize min and max coordinates with extreme values
    minCoords.x = minCoords.y = minCoords.z = PETSC_MAX_REAL;
    maxCoords.x = maxCoords.y = maxCoords.z = PETSC_MIN_REAL;
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Grid indices (Including Ghosts): xs=%d, xe=%d, ys=%d, ye=%d, zs=%d, ze=%d.\n",rank, xs, xe, ys, ye, zs, ze);
 
    // Iterate over the local grid to find min and max coordinates
    for (k = zs; k < ze; k++) {
        for (j = ys; j < ye; j++) {
            for (i = xs; i < xe; i++) {
                // Only consider nodes within the physical domain.
                if(i < user->IM && j < user->JM && k < user->KM){
                    Cmpnts coord = coordArray[k][j][i];
 
                    // Update min and max coordinates
                    if (coord.x < minCoords.x) minCoords.x = coord.x;
                    if (coord.y < minCoords.y) minCoords.y = coord.y;
                    if (coord.z < minCoords.z) minCoords.z = coord.z;
 
                    if (coord.x > maxCoords.x) maxCoords.x = coord.x;
                    if (coord.y > maxCoords.y) maxCoords.y = coord.y;
                    if (coord.z > maxCoords.z) maxCoords.z = coord.z;
                }    
            }
        }
    }
 
 
    // Add tolerance to bboxes.
    minCoords.x =  minCoords.x - BBOX_TOLERANCE;
    minCoords.y =  minCoords.y - BBOX_TOLERANCE;
    minCoords.z =  minCoords.z - BBOX_TOLERANCE;
 
    maxCoords.x =  maxCoords.x + BBOX_TOLERANCE;
    maxCoords.y =  maxCoords.y + BBOX_TOLERANCE;
    maxCoords.z =  maxCoords.z + BBOX_TOLERANCE;
 
    LOG_ALLOW(LOCAL,LOG_DEBUG," Tolerance added to the limits: %.8e .\n",(PetscReal)BBOX_TOLERANCE);
       
    // Log the computed min and max coordinates
     LOG_ALLOW(LOCAL, LOG_INFO,"[Rank %d] Bounding Box Ranges = X[%.6f, %.6f], Y[%.6f,%.6f], Z[%.6f, %.6f].\n",rank,minCoords.x, maxCoords.x,minCoords.y, maxCoords.y, minCoords.z, maxCoords.z);
 
 
    
    // Restore the coordinate array
    ierr = DMDAVecRestoreArrayRead(user->fda, coordinates, &coordArray);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error restoring coordinate array.\n");
        PROFILE_FUNCTION_END;
        return ierr;
    }
 
    // Set the local bounding box
    localBBox->min_coords = minCoords;
    localBBox->max_coords = maxCoords;
 
    // Update the bounding box inside the UserCtx for consistency
    user->bbox = *localBBox;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Exiting the function successfully.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "GatherAllBoundingBoxes"
 
/**
 * @brief Implementation of \ref GatherAllBoundingBoxes().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see GatherAllBoundingBoxes()
 */
PetscErrorCode GatherAllBoundingBoxes(UserCtx *user, BoundingBox **allBBoxes)
{
  PetscErrorCode ierr;
  PetscMPIInt    rank, size;
  BoundingBox    *bboxArray = NULL;
  BoundingBox     localBBox;
 
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  /* Validate */
  if (!user || !allBBoxes) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
                                    "GatherAllBoundingBoxes: NULL pointer");
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRMPI(ierr);
 
  /* Compute local bbox */
  ierr = ComputeLocalBoundingBox(user, &localBBox); CHKERRQ(ierr);
 
  /* Ensure everyone is synchronized before the gather */
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: about to MPI_Gather(localBBox)\n", rank);
 
  /* Allocate on root */
  if (rank == 0) {
    bboxArray = (BoundingBox*)malloc(size * sizeof(BoundingBox));
    if (!bboxArray) SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_MEM,
                                "GatherAllBoundingBoxes: malloc failed");
  }
 
  /* Collective: every rank must call */
  ierr = MPI_Gather(&localBBox, sizeof(BoundingBox), MPI_BYTE,
                    bboxArray, sizeof(BoundingBox), MPI_BYTE,
                    0, PETSC_COMM_WORLD);
  CHKERRMPI(ierr);
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: completed MPI_Gather(localBBox)\n", rank);
 
  /* Return result */
  if (rank == 0) {
    *allBBoxes = bboxArray;
  } else {
    *allBBoxes = NULL;
  }
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "BroadcastAllBoundingBoxes"
 
/**
 * @brief Internal helper implementation: `BroadcastAllBoundingBoxes()`.
 * @details Local to this translation unit.
 */
PetscErrorCode BroadcastAllBoundingBoxes(UserCtx *user, BoundingBox **bboxlist)
{
  PetscErrorCode ierr;
  (void)user;
  PetscMPIInt    rank, size;
 
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  if (!bboxlist) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
                          "BroadcastAllBoundingBoxes: NULL pointer");
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRMPI(ierr);
 
  /* Non-root ranks must allocate before the Bcast */
  if (rank != 0) {
    *bboxlist = (BoundingBox*)malloc(size * sizeof(BoundingBox));
    if (!*bboxlist) SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_MEM,
                                "BroadcastAllBoundingBoxes: malloc failed");
  }
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: about to MPI_Bcast(%d boxes)\n", rank, size);
 
  /* Collective: every rank must call */
  ierr = MPI_Bcast(*bboxlist, size * sizeof(BoundingBox), MPI_BYTE,
                   0, PETSC_COMM_WORLD);
  CHKERRMPI(ierr);
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: completed MPI_Bcast(%d boxes)\n", rank, size);
 
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "CalculateInletProperties"
/**
 * @brief Implementation of \ref CalculateInletProperties().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see CalculateInletProperties()
 */
PetscErrorCode CalculateInletProperties(UserCtx *user)
{
    PetscErrorCode ierr;
    BCFace         inlet_face_id = -1;
    PetscBool      inlet_found = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // 1. Identify the primary inlet face from the configuration
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == INLET) {
            inlet_face_id = user->boundary_faces[i].face_id;
            inlet_found = PETSC_TRUE;
            break; // Use the first inlet found
        }
    }
 
    if (!inlet_found) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "No INLET face found. Skipping inlet center calculation.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    
    Cmpnts inlet_center;
    PetscReal inlet_area;
 
    // 2. Call the generic utility to compute  the center and area of any face.
    ierr = CalculateFaceCenterAndArea(user,inlet_face_id,&inlet_center,&inlet_area); CHKERRQ(ierr);
 
    // 3. Store results in the SimCtx
    user->simCtx->CMx_c = inlet_center.x;
    user->simCtx->CMy_c = inlet_center.y;
    user->simCtx->CMz_c = inlet_center.z;
    user->simCtx->AreaInSum = inlet_area;
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Rank[%d] Inlet Center: (%.6f, %.6f, %.6f), Area: %.6f\n",
              user->simCtx->rank, inlet_center.x, inlet_center.y, inlet_center.z, inlet_area);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
 
}
 
#undef __FUNCT__
#define __FUNCT__ "CalculateOutletProperties"
/**
 * @brief Implementation of \ref CalculateOutletProperties().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see CalculateOutletProperties()
 */
PetscErrorCode CalculateOutletProperties(UserCtx *user)
{
    PetscErrorCode ierr;
    BCFace         outlet_face_id = -1;
    PetscBool      outlet_found = PETSC_FALSE;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    // 1. Identify the primary outlet face from the configuration
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == OUTLET) {
            outlet_face_id = user->boundary_faces[i].face_id;
            outlet_found = PETSC_TRUE;
            break; // Use the first outlet found
        }
    }
    if (!outlet_found) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "No OUTLET face found. Skipping outlet center calculation.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    PetscReal outlet_area;
    Cmpnts outlet_center;
    // 2. Call the generic utility to compute  the center and area of any face
    ierr = CalculateFaceCenterAndArea(user,outlet_face_id,&outlet_center,&outlet_area); CHKERRQ(ierr);
    // 3. Store results in the SimCtx
    user->simCtx->AreaOutSum = outlet_area;
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Outlet Center: (%.6f, %.6f, %.6f), Area: %.6f\n",
              outlet_center.x, outlet_center.y, outlet_center.z, outlet_area);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "CalculateFaceCenterAndArea"
/**
 * @brief Implementation of \ref CalculateFaceCenterAndArea().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/grid.h`.
 * @see CalculateFaceCenterAndArea()
 */
PetscErrorCode CalculateFaceCenterAndArea(UserCtx *user, BCFace face_id, 
                                          Cmpnts *face_center, PetscReal *face_area)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    
    // ========================================================================
    // Local accumulators for this rank's contribution
    // ========================================================================
    PetscReal      local_sum[3] = {0.0, 0.0, 0.0};  ///< Local sum of (x,y,z) coordinates
    PetscReal      localAreaSum = 0.0;               ///< Local sum of face area magnitudes
    PetscCount     local_n_points = 0;               ///< Local count of nodes
    
    // ========================================================================
    // Global accumulators after MPI reduction
    // ========================================================================
    PetscReal      global_sum[3] = {0.0, 0.0, 0.0};  ///< Global sum of coordinates
    PetscReal      globalAreaSum = 0.0;              ///< Global sum of areas
    PetscCount     global_n_points = 0;              ///< Global count of nodes
    
    // ========================================================================
    // Grid information and array pointers
    // ========================================================================
    info = user->info;
    
    // Rank's owned range in global indices
    PetscInt       xs = info.xs, xe = info.xs + info.xm;  ///< i-range: [xs, xe)
    PetscInt       ys = info.ys, ye = info.ys + info.ym;  ///< j-range: [ys, ye)
    PetscInt       zs = info.zs, ze = info.zs + info.zm;  ///< k-range: [zs, ze)
    
    // Global domain dimensions (total allocated, includes dummy at end)
    PetscInt       mx = info.mx, my = info.my, mz = info.mz;
    PetscInt       IM = user->IM;  ///< Physical domain size in i (exclude dummy)
    PetscInt       JM = user->JM;  ///< Physical domain size in j (exclude dummy)
    PetscInt       KM = user->KM;  ///< Physical domain size in k (exclude dummy)
 
    // ========================================================================
    // Interior loop bounds (adjusted to avoid ghost/boundary cells)
    // These are used for nvert checks where we need valid cell indices
    // ========================================================================
    PetscInt lxs = xs; if(xs == 0)  lxs = xs + 1;   ///< Start at 1 if on -Xi boundary
    PetscInt lxe = xe; if(xe == mx) lxe = xe - 1;   ///< End at mx-1 if on +Xi boundary
    PetscInt lys = ys; if(ys == 0)  lys = ys + 1;   ///< Start at 1 if on -Eta boundary
    PetscInt lye = ye; if(ye == my) lye = ye - 1;   ///< End at my-1 if on +Eta boundary
    PetscInt lzs = zs; if(zs == 0)  lzs = zs + 1;   ///< Start at 1 if on -Zeta boundary
    PetscInt lze = ze; if(ze == mz) lze = ze - 1;   ///< End at mz-1 if on +Zeta boundary
 
    // ========================================================================
    // Physical node bounds (exclude dummy indices at mx-1, my-1, mz-1)
    // These are used for coordinate loops where we want ALL physical nodes
    // ========================================================================
    PetscInt i_max = (xe == mx) ? mx - 1 : xe;  ///< Exclude dummy at i=mx-1 (e.g., i=25)
    PetscInt j_max = (ye == my) ? my - 1 : ye;  ///< Exclude dummy at j=my-1 (e.g., j=25)
    PetscInt k_max = (ze == mz) ? mz - 1 : ze;  ///< Exclude dummy at k=mz-1 (e.g., k=97)
 
    // ========================================================================
    // Array pointers for field access
    // ========================================================================
    Vec            lCoor;                       ///< Local ghosted coordinate vector
    Cmpnts         ***coor;                    ///< Nodal coordinates [k][j][i]
    Cmpnts         ***csi, ***eta, ***zet;    ///< Face metric tensors [k][j][i]
    PetscReal      ***nvert;                   ///< Cell blanking field [k][j][i] (shifted +1)
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // ========================================================================
    // Step 1: Check if this rank owns the specified boundary face
    // ========================================================================
    PetscBool owns_face = PETSC_FALSE;
    ierr = CanRankServiceFace(&info,IM,JM,KM,face_id,&owns_face); CHKERRQ(ierr);
    if(owns_face){    
        // ========================================================================
        // Step 2: Get read-only array access for all required fields
        // ========================================================================
        ierr = DMGetCoordinatesLocal(user->da, &lCoor); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, lCoor, &coor); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->da, user->lNvert, &nvert); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
 
        // ========================================================================
        // Step 3: Loop over the specified face and accumulate center and area
        // ========================================================================
        switch (face_id) {
            
            // ====================================================================
            // BC_FACE_NEG_X: Face at i=0 (bottom boundary in i-direction)
            // ====================================================================
            case BC_FACE_NEG_X:
                if (xs == 0) {
                    PetscInt i = 0;  // Face is at node index i=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // Loop over ALL physical nodes on this face
                    // For my=26, mz=98: j∈[0,24], k∈[0,96] → 25×97 = 2,425 nodes
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt j = ys; j < j_max; j++) {
                            // Accumulate coordinates at node [k][j][0]
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // Loop over interior range where nvert checks are valid
                    // For my=26, mz=98: j∈[1,24], k∈[1,96] → 24×96 = 2,304 cells
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt j = lys; j < lye; j++) {
                            // Check if adjacent cell is fluid
                            // nvert[k][j][i+1] = nvert[k][j][1] checks Cell 0
                            // (Physical Cell 0 in j-k plane, stored at shifted index [1])
                            if (nvert[k][j][i+1] < 0.1) {
                                // Cell is fluid - add face area contribution
                                // Face area = magnitude of csi metric at [k][j][0]
                                localAreaSum += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                                csi[k][j][i].y * csi[k][j][i].y +
                                                csi[k][j][i].z * csi[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_X: Face at i=IM-1 (top boundary in i-direction)
            // ====================================================================
            case BC_FACE_POS_X:
                if (xe == mx) {
                    PetscInt i = mx - 2;  // Last physical node (e.g., i=24 for mx=26)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt j = ys; j < j_max; j++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt j = lys; j < lye; j++) {
                            // Check if adjacent cell is fluid
                            // nvert[k][j][i] = nvert[k][j][24] checks last cell (Cell 23)
                            // (Physical Cell 23, stored at shifted index [24])
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of csi metric at [k][j][24]
                                localAreaSum += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                                csi[k][j][i].y * csi[k][j][i].y +
                                                csi[k][j][i].z * csi[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_NEG_Y: Face at j=0 (bottom boundary in j-direction)
            // ====================================================================
            case BC_FACE_NEG_Y:
                if (ys == 0) {
                    PetscInt j = 0;  // Face is at node index j=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, mz=98: i∈[0,24], k∈[0,96] → 25×97 = 2,425 nodes
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, mz=98: i∈[1,24], k∈[1,96] → 24×96 = 2,304 cells
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j+1][i] = nvert[k][1][i] checks Cell 0
                            if (nvert[k][j+1][i] < 0.1) {
                                // Face area = magnitude of eta metric at [k][0][i]
                                localAreaSum += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                                eta[k][j][i].y * eta[k][j][i].y +
                                                eta[k][j][i].z * eta[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_Y: Face at j=JM-1 (top boundary in j-direction)
            // ====================================================================
            case BC_FACE_POS_Y:
                if (ye == my) {
                    PetscInt j = my - 2;  // Last physical node (e.g., j=24 for my=26)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j][i] = nvert[k][24][i] checks last cell (Cell 23)
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of eta metric at [k][24][i]
                                localAreaSum += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                                eta[k][j][i].y * eta[k][j][i].y +
                                                eta[k][j][i].z * eta[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_NEG_Z: Face at k=0 (inlet, bottom boundary in k-direction)
            // ====================================================================
            case BC_FACE_NEG_Z:
                if (zs == 0) {
                    PetscInt k = 0;  // Face is at node index k=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, my=26: i∈[0,24], j∈[0,24] → 25×25 = 625 nodes
                    for (PetscInt j = ys; j < j_max; j++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, my=26: i∈[1,24], j∈[1,24] → 24×24 = 576 cells
                    for (PetscInt j = lys; j < lye; j++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k+1][j][i] = nvert[1][j][i] checks Cell 0
                            // (Physical Cell 0 in i-j plane, stored at shifted index [1])
                            if (nvert[k+1][j][i] < 0.1) {
                                // Face area = magnitude of zet metric at [0][j][i]
                                localAreaSum += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                                zet[k][j][i].y * zet[k][j][i].y +
                                                zet[k][j][i].z * zet[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_Z: Face at k=KM-1 (outlet, top boundary in k-direction)
            // ====================================================================
            case BC_FACE_POS_Z:
                if (ze == mz) {
                    PetscInt k = mz - 2;  // Last physical node (e.g., k=96 for mz=98)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, my=26: i∈[0,24], j∈[0,24] → 25×25 = 625 nodes
                    for (PetscInt j = ys; j < j_max; j++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, my=26: i∈[1,24], j∈[1,24] → 24×24 = 576 cells
                    for (PetscInt j = lys; j < lye; j++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j][i] = nvert[96][j][i] checks last cell (Cell 95)
                            // (Physical Cell 95, stored at shifted index [96])
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of zet metric at [96][j][i]
                                localAreaSum += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                                zet[k][j][i].y * zet[k][j][i].y +
                                                zet[k][j][i].z * zet[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            default:
                SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 
                        "Unknown face_id %d in CalculateFaceCenterAndArea", face_id);
        }
        
        // ========================================================================
        // Step 4: Restore array access (release pointers)
        // ========================================================================
        ierr = DMDAVecRestoreArrayRead(user->fda, lCoor, &coor); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, &nvert); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    }    
    // ========================================================================
    // Step 5: Perform MPI reductions to get global sums
    // ========================================================================
    // Sum coordinate contributions from all ranks
    ierr = MPI_Allreduce(local_sum, global_sum, 3, MPI_DOUBLE, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // Sum node counts from all ranks
    ierr = MPI_Allreduce(&local_n_points, &global_n_points, 1, MPI_COUNT, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // Sum area contributions from all ranks
    ierr = MPI_Allreduce(&localAreaSum, &globalAreaSum, 1, MPI_DOUBLE, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    // ========================================================================
    // Step 6: Calculate geometric center by averaging coordinates
    // ========================================================================
    if (global_n_points > 0) {
        face_center->x = global_sum[0] / global_n_points;
        face_center->y = global_sum[1] / global_n_points;
        face_center->z = global_sum[2] / global_n_points;
        LOG_ALLOW(GLOBAL, LOG_INFO, 
                  "Calculated center for Face %s: (x=%.4f, y=%.4f, z=%.4f) from %lld nodes\n", 
                  BCFaceToString(face_id), 
                  face_center->x, face_center->y, face_center->z,
                  (long long)global_n_points);
    } else {
        // No nodes found - this should not happen for a valid face
        LOG_ALLOW(GLOBAL, LOG_WARNING, 
                  "WARNING: Face %s identified but no grid points found. Center not calculated.\n",
                  BCFaceToString(face_id));
        face_center->x = face_center->y = face_center->z = 0.0;
    }
    
    // ========================================================================
    // Step 7: Return computed total area
    // ========================================================================
    *face_area = globalAreaSum;
    LOG_ALLOW(GLOBAL, LOG_INFO, 
              "Calculated area for Face %s: Area=%.6f\n", 
              BCFaceToString(face_id), *face_area);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}