interpolation.c

Go to the documentation of this file.

/**
 * @file interpolation.c
 * @brief Main program for DMSwarm interpolation using the fdf-curvIB method.
 *
 * Provides routines for interpolation between corner-based and center-based
 * fields in the cell-centered DM (fda), plus partial usage examples for
 * DMSwarm-based field sampling.
 */
 
#include "interpolation.h"
 
// Number of weights used in certain trilinear interpolation examples
#define NUM_WEIGHTS 8
// Define a buffer size for error messages if not already available
#ifndef ERROR_MSG_BUFFER_SIZE
#define ERROR_MSG_BUFFER_SIZE 256 // Or use PETSC_MAX_PATH_LEN if appropriate
#endif
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateFieldFromCornerToCenter_Vector"
/**
 * @brief Internal helper implementation: `InterpolateFieldFromCornerToCenter_Vector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateFieldFromCornerToCenter_Vector(
    Cmpnts ***field_arr,     /* Input: Ghosted local array view from user->fda (global node indices) */
    Cmpnts ***centfield_arr, /* Output: Array view for cell-centered data (global indices) */
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i].x = 0.125 * (field_arr[k][j][i].x + field_arr[k][j-1][i].x +
                                                    field_arr[k-1][j][i].x + field_arr[k-1][j-1][i].x +
                                                    field_arr[k][j][i-1].x + field_arr[k][j-1][i-1].x +
                                                    field_arr[k-1][j][i-1].x + field_arr[k-1][j-1][i-1].x);
 
                centfield_arr[k][j][i].y = 0.125 * (field_arr[k][j][i].y + field_arr[k][j-1][i].y +
                                                    field_arr[k-1][j][i].y + field_arr[k-1][j-1][i].y +
                                                    field_arr[k][j][i-1].y + field_arr[k][j-1][i-1].y +
                                                    field_arr[k-1][j][i-1].y + field_arr[k-1][j-1][i-1].y);
 
                centfield_arr[k][j][i].z = 0.125 * (field_arr[k][j][i].z + field_arr[k][j-1][i].z +
                                                    field_arr[k-1][j][i].z + field_arr[k-1][j-1][i].z +
                                                    field_arr[k][j][i-1].z + field_arr[k][j-1][i-1].z +
                                                    field_arr[k-1][j][i-1].z + field_arr[k-1][j-1][i-1].z);
            }
        }
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateFieldFromCornerToCenter_Scalar"
/**
 * @brief Internal helper implementation: `InterpolateFieldFromCornerToCenter_Scalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateFieldFromCornerToCenter_Scalar(
    PetscReal ***field_arr,     /* Input: Ghosted local array view from user->fda (global node indices) */
    PetscReal ***centfield_arr, /* Output: Array view for cell-centered data (global indices) */
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i] = 0.125 * (field_arr[k][j][i] + field_arr[k][j-1][i] +
                                                  field_arr[k-1][j][i] + field_arr[k-1][j-1][i] +
                                                  field_arr[k][j][i-1] + field_arr[k][j-1][i-1] +
                                                  field_arr[k-1][j][i-1] + field_arr[k-1][j-1][i-1]);
            }
        }
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "TestCornerToCenterInterpolation"
/**
 * @brief Internal helper implementation: `TestCornerToCenterInterpolation()`.
 * @details Local to this translation unit.
 */
PetscErrorCode TestCornerToCenterInterpolation(UserCtx *user)
{
    PetscErrorCode ierr;
    Vec            lCoords, TestCent;
    Cmpnts         ***coor_arr, ***test_cent_arr;
    PetscReal      diff_norm;
 
    PetscFunctionBeginUser;
 
    // 1. Create a temporary vector to hold the result of our interpolation.
    //    It must have the same layout and size as the ground-truth user->Cent vector.
    ierr = VecDuplicate(user->Cent, &TestCent); CHKERRQ(ierr);
 
    // 2. Get the input (corner coordinates) and output (our test vector) arrays.
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 3. Call the generic interpolation macro.
    //    The macro will see that `coor_arr` is of type `Cmpnts***` and correctly
    //    call the `InterpolateFieldFromCornerToCenter_Vector` function.
    ierr = InterpolateFieldFromCornerToCenter(coor_arr, test_cent_arr, user); CHKERRQ(ierr);
 
    // 4. Restore the arrays.
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 5. IMPORTANT: Assemble the vector so its values are communicated across processors
    //    and it's ready for global operations like VecNorm.
    ierr = VecAssemblyBegin(TestCent); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(TestCent); CHKERRQ(ierr);
 
    // 6. Compare the result with the ground truth.
    //    We compute TestCent = -1.0 * user->Cent + 1.0 * TestCent.
    //    This calculates the difference vector: TestCent - user->Cent.
    ierr = VecAXPY(TestCent, -1.0, user->Cent); CHKERRQ(ierr);
 
    //    Now, compute the L2 norm of the difference vector. If the functions are
    //    identical, the norm should be zero (or very close due to floating point).
    ierr = VecNorm(TestCent, NORM_2, &diff_norm); CHKERRQ(ierr);
 
    // 7. Report the result and clean up.
    if (diff_norm < 1.0e-12) {
        LOG_ALLOW(GLOBAL,LOG_DEBUG,"[SUCCESS] Test passed. Norm of difference is %g.\n", (double)diff_norm);
    } else {
        LOG_ALLOW(GLOBAL,LOG_DEBUG, "[FAILURE] Test failed. Norm of difference is %g.\n", (double)diff_norm);
    }
 
    ierr = VecDestroy(&TestCent); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateFieldFromCenterToCorner_Vector"
/**
 * @brief Internal helper implementation: `InterpolateFieldFromCenterToCorner_Vector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateFieldFromCenterToCorner_Vector(
    Cmpnts ***centfield_arr, /* Input: Ghosted local array from Vec (read) */
    Cmpnts ***corner_arr,    /* Output: global array from Vec (write) */
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                Cmpnts sum = {0.0, 0.0, 0.0};
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
                    
                                    sum.x += cell_val.x;
                                    sum.y += cell_val.y;
                                    sum.z += cell_val.z;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateFieldFromCenterToCorner_Scalar"
/**
 * @brief Internal helper implementation: `InterpolateFieldFromCenterToCorner_Scalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar(
    PetscReal ***centfield_arr, /* Input: Ghosted local array from Vec (read) */
    PetscReal ***corner_arr,    /* Output: global array from Vec (write) */
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                PetscReal sum = 0.0;
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    PetscReal cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val);
                    
                                    sum += cell_val;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write] = sum / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = 0.0;
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}
 
#undef __FUNCT
#define __FUNCT "PiecWiseLinearInterpolation_Scalar"
/**
 * @brief Internal helper implementation: `PieceWiseLinearInterpolation_Scalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode PieceWiseLinearInterpolation_Scalar(
    const char   *fieldName,
    PetscReal  ***fieldScal,
    PetscInt      iCell,
    PetscInt      jCell,
    PetscInt      kCell,
    PetscReal    *val)
{
  PetscFunctionBegin;
  *val = fieldScal[kCell][jCell][iCell];
  
  // Optional logging
  LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "Field '%s' at (i=%d, j=%d, k=%d) => val=%.6f\n",
      fieldName, iCell, jCell, kCell, *val);
 
  PetscFunctionReturn(0);
}
 
#undef __FUNCT
#define __FUNCT "PiecWiseLinearInterpolation_Vector"
/**
 * @brief Internal helper implementation: `PieceWiseLinearInterpolation_Vector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode PieceWiseLinearInterpolation_Vector(
    const char   *fieldName,
    Cmpnts     ***fieldVec,
    PetscInt      iCell,
    PetscInt      jCell,
    PetscInt      kCell,
    Cmpnts       *vec)
{
  PetscFunctionBegin;
  vec->x = fieldVec[kCell][jCell][iCell].x;
  vec->y = fieldVec[kCell][jCell][iCell].y;
  vec->z = fieldVec[kCell][jCell][iCell].z;
 
  // Optional logging
  LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "Field '%s' at (i=%d, j=%d, k=%d) => (x=%.6f, y=%.6f, z=%.6f)\n",
      fieldName, iCell, jCell, kCell, vec->x, vec->y, vec->z);
 
  PetscFunctionReturn(0);
}
 
 
#undef __FUNCT
#define __FUNCT "ComputeTrilinearWeights"
 
/**
 * @brief Internal helper implementation: `ComputeTrilinearWeights()`.
 * @details Local to this translation unit.
 */
static inline void ComputeTrilinearWeights(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w) {
    LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
 
    // Ensure a1, a2, a3 are within [0,1]
    a1 = PetscMax(0.0, PetscMin(1.0, a1));
    a2 = PetscMax(0.0, PetscMin(1.0, a2));
    a3 = PetscMax(0.0, PetscMin(1.0, a3));
 
    const PetscReal oa1 = 1.0 - a1;
    const PetscReal oa2 = 1.0 - a2;
    const PetscReal oa3 = 1.0 - a3;
 
    w[0] = oa1 * oa2 * oa3;  /* cornerOffsets[0] => (0,0,0) */
    w[1] = a1  * oa2 * oa3;  /* cornerOffsets[1] => (1,0,0) */
    w[2] = oa1 * a2  * oa3;  /* cornerOffsets[2] => (0,1,0) */
    w[3] = a1  * a2  * oa3;  /* cornerOffsets[3] => (1,1,0) */
    w[4] = oa1 * oa2 * a3;   /* cornerOffsets[4] => (0,0,1) */
    w[5] = a1  * oa2 * a3;   /* cornerOffsets[5] => (1,0,1) */
    w[6] = oa1 * a2  * a3;   /* cornerOffsets[6] => (0,1,1) */
    w[7] = a1  * a2  * a3;   /* cornerOffsets[7] => (1,1,1) */
 
    // Log the computed weights for debugging
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "Weights computed - "
        "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f. \n",
        w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
}
 
/**
 * @brief Unclamped trilinear weights for boundary extrapolation.
 * @details Identical to ComputeTrilinearWeights() but without clamping
 *          a1, a2, a3 to [0,1]. Allows weights outside [0,1] for linear
 *          extrapolation at non-periodic boundaries. Weights still sum to 1.0.
 */
static inline void ComputeTrilinearWeightsUnclamped(PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *w) {
    LOG_ALLOW(GLOBAL, LOG_VERBOSE, "Computing unclamped weights for a1=%f, a2=%f, a3=%f.\n", a1, a2, a3);
 
    const PetscReal oa1 = 1.0 - a1;
    const PetscReal oa2 = 1.0 - a2;
    const PetscReal oa3 = 1.0 - a3;
 
    w[0] = oa1 * oa2 * oa3;
    w[1] = a1  * oa2 * oa3;
    w[2] = oa1 * a2  * oa3;
    w[3] = a1  * a2  * oa3;
    w[4] = oa1 * oa2 * a3;
    w[5] = a1  * oa2 * a3;
    w[6] = oa1 * a2  * a3;
    w[7] = a1  * a2  * a3;
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE, "Unclamped weights - "
        "w0=%f, w1=%f, w2=%f, w3=%f, w4=%f, w5=%f, w6=%f, w7=%f.\n",
        w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7]);
}
 
#undef __FUNCT
#define __FUNCT "TrilinearInterpolation_Scalar"
 
/**
 * @brief Internal helper implementation: `TrilinearInterpolation_Scalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode TrilinearInterpolation_Scalar(
    const char   *fieldName,
    PetscReal  ***fieldScal,
    PetscInt      i,
    PetscInt      j,
    PetscInt      k,
    PetscReal     a1,
    PetscReal     a2,
    PetscReal     a3,
    PetscReal    *val)
{
    PetscFunctionBegin; // PETSc macro for error/stack tracing
 
    // Compute the 8 corner weights
    PetscReal wcorner[8];
    ComputeTrilinearWeights(a1, a2, a3, wcorner);
 
    // Offsets for cell corners
    PetscInt i1 = i + 1;
    PetscInt j1 = j + 1;
    PetscInt k1 = k + 1;
 
    // Initialize the output scalar
    PetscReal sum = 0.0;
 
    // Corner 0 => (i, j, k)
    sum += wcorner[0] * fieldScal[k ][j ][i ];
    // Corner 1 => (i+1, j, k)
    sum += wcorner[1] * fieldScal[k ][j ][i1];
    // Corner 2 => (i, j+1, k)
    sum += wcorner[2] * fieldScal[k ][j1][i ];
    // Corner 3 => (i+1, j+1, k)
    sum += wcorner[3] * fieldScal[k ][j1][i1];
    // Corner 4 => (i, j, k+1)
    sum += wcorner[4] * fieldScal[k1][j ][i ];
    // Corner 5 => (i+1, j, k+1)
    sum += wcorner[5] * fieldScal[k1][j ][i1];
    // Corner 6 => (i, j+1, k+1)
    sum += wcorner[6] * fieldScal[k1][j1][i ];
    // Corner 7 => (i+1, j+1, k+1)
    sum += wcorner[7] * fieldScal[k1][j1][i1];
 
    *val = sum;
 
    // Logging (optional)
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
        "Field '%s' at (i=%d, j=%d, k=%d), "
        "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
        fieldName, i, j, k, a1, a2, a3, *val);
 
    // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
    //    "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
      //    fieldName);
 
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT
#define __FUNCT "TrilinearInterpolation_Vector"
/**
 * @brief Internal helper implementation: `TrilinearInterpolation_Vector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode TrilinearInterpolation_Vector(
    const char   *fieldName,
    Cmpnts     ***fieldVec,  /* 3D array [k][j][i], dimension [mz][my][mx] */
    PetscInt      i, // local cell index i
    PetscInt      j, // local cell index j
    PetscInt      k, // local cell index k
    PetscReal     a1,
    PetscReal     a2,
    PetscReal     a3,
    Cmpnts       *vec)
{
  PetscFunctionBegin; // PETSc macro for error/stack tracing
 
  // Compute the 8 corner weights
  PetscErrorCode ierr;
  PetscReal wcorner[8];
  PetscMPIInt rank;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
  ComputeTrilinearWeights(a1, a2, a3, wcorner);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d]  Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
 
  // For partial interpolation, we'll keep track of how many corners are valid
  // and how much sum of weights is used. Then we do a final normalization.
  PetscReal sumW = 0.0;
  Cmpnts    accum = {0.0, 0.0, 0.0};
 
  // The eight corner indices, with their weights:
  // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
  // We store them in an array to iterate cleanly.
  const PetscInt cornerOffsets[8][3] = {
    {0, 0, 0},
    {1, 0, 0},
    {0, 1, 0},
    {1, 1, 0},
    {0, 0, 1},
    {1, 0, 1},
    {0, 1, 1},
    {1, 1, 1}
  };
 
  // Weighted partial sum
  for (PetscInt c = 0; c < 8; c++) {
    const PetscInt di = cornerOffsets[c][0];
    const PetscInt dj = cornerOffsets[c][1];
    const PetscInt dk = cornerOffsets[c][2];
    PetscInt iC = i + di;
    PetscInt jC = j + dj;
    PetscInt kC = k + dk;
 
    /*
    // skip if out of domain
    // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
    if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
        jC < 0 || jC >= (PetscInt)userGlobalMy ||
        kC < 0 || kC >= (PetscInt)userGlobalMz)
    {
      // skip this corner
      continue;
    }
 
    */
 
    LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] %s[%d][%d][%d] = (%.4f,%.4f,%.4f).\n",rank,fieldName,kC,jC,iC,fieldVec[kC][jC][iC].x,fieldVec[kC][jC][iC].y,fieldVec[kC][jC][iC].z);
    
    // Otherwise, accumulate
    accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
    accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
    accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
    sumW += wcorner[c];
  }
 
  // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
  if (sumW > 1.0e-14) {
    vec->x = accum.x / sumW;
    vec->y = accum.y / sumW;
    vec->z = accum.z / sumW;
  } else {
    vec->x = 0.0;  vec->y = 0.0;  vec->z = 0.0;
  }
 
  PetscFunctionReturn(0);
}
 
 
#undef __FUNCT
#define __FUNCT "InterpolateEulerFieldToSwarmForParticle"
/**
 * @brief Internal helper implementation: `InterpolateEulerFieldToSwarmForParticle()`.
 * @details Local to this translation unit.
 */
static inline PetscErrorCode InterpolateEulerFieldToSwarmForParticle(
    const char  *fieldName,
    void        *fieldPtr,   /* typed Pointer => either (PetscReal***) or (Cmpnts***) */
    Particle    *particle,   /* particle struct containing Cell ID and weight information. */
    void        *swarmOut,   /* typed Pointer => (PetscReal*) or (Cmpnts*) or dof=3 array */
    PetscInt     p,          /* particle index */
    PetscInt     blockSize)  /* dof=1 => scalar, dof=3 => vector */
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
 
  PetscInt iCell = particle->cell[0];
  PetscInt jCell = particle->cell[1];
  PetscInt kCell = particle->cell[2];
  PetscReal a1 = particle->weights.x;
  PetscReal a2 = particle->weights.y;
  PetscReal a3 = particle->weights.z;
 
  PROFILE_FUNCTION_BEGIN;
 
  // Optional logging at start
  LOG_ALLOW(LOCAL, LOG_VERBOSE,
    "field='%s', blockSize=%d, "
    "cell IDs=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
    fieldName, blockSize, iCell, jCell, kCell, a1, a2, a3);
 
  /*
     If blockSize=1, we PetscInterpret the fieldPtr as a 3D array of PetscReal (scalar).
     If blockSize=3, we PetscInterpret the fieldPtr as a 3D array of Cmpnts (vector).
   */
  if (blockSize == 1) {
    /* Scalar field: Cast fieldPtr to (PetscReal ***). */
    PetscReal ***fieldScal = (PetscReal ***) fieldPtr;
    PetscReal val;
 
    // Currently using trilinear.
    ierr = TrilinearInterpolation(fieldName, fieldScal,
                                   iCell, jCell, kCell,
                                   a1, a2, a3,
                                   &val);
     CHKERRQ(ierr);
 
    // Alternative (commented) call to PiecewiseLinearInterpolation (zeroth order) :
    //   PetscErrorCode ierr = PiecewiseLinearInterpolation(fieldName,
    //                                                   fieldScal,
    //                                                   iCell, jCell, kCell,
    //                                                   &val);
    // CHKERRQ(ierr);
 
    // Write the scalar result to the swarm output at index [p].
    ((PetscReal*)swarmOut)[p] = val;
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "field='%s', result=%.6f "
      "stored at swarmOut index p=%d.\n", fieldName, val, (PetscInt)p);
  }
  else if (blockSize == 3) {
    /* Vector field: Cast fieldPtr to (Cmpnts ***). */
    Cmpnts ***fieldVec = (Cmpnts ***) fieldPtr;
    Cmpnts vec;
 
 
    
    
    // Piecewise interpolation (zeroth order).
    //  PetscErrorCode ierr = PieceWiseLinearInterpolation(fieldName,
    //                                                   fieldVec,
    //                                                   iCell, jCell, kCell,
    //                                                   &vec);
    // CHKERRQ(ierr);
 
    // Alternative (commented) call to trilinear:
     ierr = TrilinearInterpolation(fieldName, fieldVec,
                                   iCell, jCell, kCell,
                                   a1, a2, a3,
                                   &vec);
     CHKERRQ(ierr);
 
    // If swarmOut is an array of 3 reals per particle:
    ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
    ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
    ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
      "field='%s', result=(%.6f,%.6f,%.6f) "
      "stored at swarmOut[3p..3p+2], p=%d.\n",
      fieldName, vec.x, vec.y, vec.z, (PetscInt)p);
 
    /*
       If you store the vector result as a Cmpnts in the swarm, do instead:
         ((Cmpnts*)swarmOut)[p] = vec;
       but ensure your DMSwarm field is sized for a Cmpnts struct.
    */
  }
  else {
    /* If blockSize isn't 1 or 3, we raise an error. */
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
      "InterpolateEulerFieldToSwarmForParticle: only blockSize=1 or 3 supported, got %d.",
      (PetscInt)blockSize);
  }
 
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateEulerFieldFromCenterToSwarm"
 
/**
 * @brief Direct cell-center trilinear interpolation (second-order on curvilinear grids).
 * @details For each particle, determines the 8 nearest cell centers via octant
 *          detection, constructs a dual cell, computes face-distance-based
 *          trilinear weights, and interpolates directly from cell-centered data.
 *          At non-periodic boundaries where the dual cell cannot be fully formed,
 *          octant clamping with unclamped trilinear extrapolation preserves
 *          second-order accuracy. Periodic boundaries use ghost cell data directly.
 *          No intermediate corner staging or extra ghost exchange is needed.
 */
static PetscErrorCode InterpolateEulerFieldFromCenterToSwarm(
    UserCtx    *user,
    Vec         fieldLocal_cellCentered,
    const char *fieldName,
    const char *swarmOutFieldName)
{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    PetscInt       bs;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    void          *fieldPtr = NULL;
    Cmpnts       ***cent = NULL;
 
    PetscInt  *cellIDs = NULL;
    PetscReal *weights = NULL;
    PetscInt64 *pids = NULL;
    void      *swarmOut = NULL;
    PetscReal *pos = NULL;
    PetscInt  *status = NULL;
    PetscInt   nLocal;
 
    PetscFunctionBegin;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
    if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting direct cell-center interpolation for field '%s'...\n", fieldName);
 
    /* Number of physical cells in each direction */
    PetscInt nCellsX = info.mx - 2;
    PetscInt nCellsY = info.my - 2;
    PetscInt nCellsZ = info.mz - 2;
 
    /* Periodic flags per direction */
    PetscBool x_periodic = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
    PetscBool y_periodic = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
    PetscBool z_periodic = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
 
    /* Get read-only arrays for the cell-centered field and cell center coordinates */
    DM dm_field = (bs == 3) ? user->fda : user->da;
    ierr = DMDAVecGetArrayRead(dm_field, fieldLocal_cellCentered, &fieldPtr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCent, (void *)&cent); CHKERRQ(ierr);
 
    /* Retrieve swarm fields */
    ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
 
    /* Loop over each local particle */
    for (PetscInt p = 0; p < nLocal; p++) {
 
        Particle particle;
        ierr = UnpackSwarmFields(p, pids, weights, pos, cellIDs, NULL, status, NULL, NULL, NULL, &particle); CHKERRQ(ierr);
 
        /* Step 1: Get host cell indices and existing weights */
        PetscInt ci = particle.cell[0];
        PetscInt cj = particle.cell[1];
        PetscInt ck = particle.cell[2];
        PetscReal a1 = particle.weights.x;
        PetscReal a2 = particle.weights.y;
        PetscReal a3 = particle.weights.z;
 
        /* Step 2: Determine octant */
        PetscInt oi = (a1 < 0.5) ? -1 : 0;
        PetscInt oj = (a2 < 0.5) ? -1 : 0;
        PetscInt ok = (a3 < 0.5) ? -1 : 0;
        PetscInt bi = ci + oi;
        PetscInt bj = cj + oj;
        PetscInt bk = ck + ok;
 
        /* Step 3: Clamp octant at non-periodic boundaries */
        PetscBool needs_unclamped = PETSC_FALSE;
 
        if (!x_periodic) {
            if (bi < 0)              { bi = 0;            needs_unclamped = PETSC_TRUE; }
            if (bi + 1 >= nCellsX)   { bi = nCellsX - 2;  needs_unclamped = PETSC_TRUE; }
        }
        if (!y_periodic) {
            if (bj < 0)              { bj = 0;            needs_unclamped = PETSC_TRUE; }
            if (bj + 1 >= nCellsY)   { bj = nCellsY - 2;  needs_unclamped = PETSC_TRUE; }
        }
        if (!z_periodic) {
            if (bk < 0)              { bk = 0;            needs_unclamped = PETSC_TRUE; }
            if (bk + 1 >= nCellsZ)   { bk = nCellsZ - 2;  needs_unclamped = PETSC_TRUE; }
        }
 
        /* Step 4: Bounds check — dual cell stencil must fit in ghosted region */
        if (bi + 1 < info.gxs || bi + 1 >= info.gxs + info.gxm - 1 ||
            bj + 1 < info.gys || bj + 1 >= info.gys + info.gym - 1 ||
            bk + 1 < info.gzs || bk + 1 >= info.gzs + info.gzm - 1)
        {
            LOG_ALLOW(LOCAL, LOG_WARNING,
                "[Rank %d] Particle PID %lld: dual cell (%d,%d,%d) out of ghosted region. Zeroing '%s'.\n",
                rank, (long long)particle.PID, bi, bj, bk, fieldName);
            if (bs == 3) {
                ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
            } else {
                ((PetscReal*)swarmOut)[p] = 0.0;
            }
            continue;
        }
 
        /* Step 5: Construct dual cell from lCent (shifted index: cent[k+1][j+1][i+1] for cell (i,j,k))
         * Vertex ordering follows GetCellVerticesFromGrid() convention (walkingsearch.c:423-430). */
        Cell dual_cell;
        dual_cell.vertices[0] = cent[bk + 1][bj + 1][bi + 1];   /* (bi,   bj,   bk  ) */
        dual_cell.vertices[1] = cent[bk + 1][bj + 1][bi + 2];   /* (bi+1, bj,   bk  ) */
        dual_cell.vertices[2] = cent[bk + 1][bj + 2][bi + 2];   /* (bi+1, bj+1, bk  ) */
        dual_cell.vertices[3] = cent[bk + 1][bj + 2][bi + 1];   /* (bi,   bj+1, bk  ) */
        dual_cell.vertices[4] = cent[bk + 2][bj + 2][bi + 1];   /* (bi,   bj+1, bk+1) */
        dual_cell.vertices[5] = cent[bk + 2][bj + 2][bi + 2];   /* (bi+1, bj+1, bk+1) */
        dual_cell.vertices[6] = cent[bk + 2][bj + 1][bi + 2];   /* (bi+1, bj,   bk+1) */
        dual_cell.vertices[7] = cent[bk + 2][bj + 1][bi + 1];   /* (bi,   bj,   bk+1) */
 
        /* Step 6: Compute face distances from particle to dual cell faces */
        PetscReal d[NUM_FACES];
        ierr = CalculateDistancesToCellFaces(particle.loc, &dual_cell, d, 1e-11); CHKERRQ(ierr);
 
        /* Step 7: Compute trilinear parametric coordinates from face distances */
        PetscReal a1_new = d[LEFT]   / (d[LEFT]   + d[RIGHT]);
        PetscReal a2_new = d[BOTTOM] / (d[BOTTOM] + d[TOP]);
        PetscReal a3_new = d[BACK]   / (d[FRONT]  + d[BACK]);
 
        LOG_ALLOW(LOCAL, LOG_VERBOSE,
            "[Rank %d] PID %lld: dual base=(%d,%d,%d), new weights=(%.4f,%.4f,%.4f), unclamped=%d\n",
            rank, (long long)particle.PID, bi, bj, bk, a1_new, a2_new, a3_new, (int)needs_unclamped);
 
        /* Step 8: Call TrilinearInterpolation directly with shifted dual-cell indices.
         * Shifted index: pass (bi+1, bj+1, bk+1) so the kernel accesses
         * field[bk+1..bk+2][bj+1..bj+2][bi+1..bi+2] = cells (bi..bi+1, bj..bj+1, bk..bk+1). */
        if (bs == 1) {
            PetscReal ***fieldScal = (PetscReal ***)fieldPtr;
            PetscReal val;
            if (needs_unclamped) {
                PetscReal w[8];
                ComputeTrilinearWeightsUnclamped(a1_new, a2_new, a3_new, w);
                val = w[0] * fieldScal[bk+1][bj+1][bi+1] + w[1] * fieldScal[bk+1][bj+1][bi+2]
                    + w[2] * fieldScal[bk+1][bj+2][bi+1] + w[3] * fieldScal[bk+1][bj+2][bi+2]
                    + w[4] * fieldScal[bk+2][bj+1][bi+1] + w[5] * fieldScal[bk+2][bj+1][bi+2]
                    + w[6] * fieldScal[bk+2][bj+2][bi+1] + w[7] * fieldScal[bk+2][bj+2][bi+2];
            } else {
                ierr = TrilinearInterpolation_Scalar(fieldName, fieldScal,
                    bi + 1, bj + 1, bk + 1, a1_new, a2_new, a3_new, &val); CHKERRQ(ierr);
            }
            ((PetscReal*)swarmOut)[p] = val;
        } else {
            Cmpnts ***fieldVec = (Cmpnts ***)fieldPtr;
            Cmpnts vec;
            if (needs_unclamped) {
                PetscReal w[8];
                ComputeTrilinearWeightsUnclamped(a1_new, a2_new, a3_new, w);
                vec.x = w[0]*fieldVec[bk+1][bj+1][bi+1].x + w[1]*fieldVec[bk+1][bj+1][bi+2].x
                      + w[2]*fieldVec[bk+1][bj+2][bi+1].x + w[3]*fieldVec[bk+1][bj+2][bi+2].x
                      + w[4]*fieldVec[bk+2][bj+1][bi+1].x + w[5]*fieldVec[bk+2][bj+1][bi+2].x
                      + w[6]*fieldVec[bk+2][bj+2][bi+1].x + w[7]*fieldVec[bk+2][bj+2][bi+2].x;
                vec.y = w[0]*fieldVec[bk+1][bj+1][bi+1].y + w[1]*fieldVec[bk+1][bj+1][bi+2].y
                      + w[2]*fieldVec[bk+1][bj+2][bi+1].y + w[3]*fieldVec[bk+1][bj+2][bi+2].y
                      + w[4]*fieldVec[bk+2][bj+1][bi+1].y + w[5]*fieldVec[bk+2][bj+1][bi+2].y
                      + w[6]*fieldVec[bk+2][bj+2][bi+1].y + w[7]*fieldVec[bk+2][bj+2][bi+2].y;
                vec.z = w[0]*fieldVec[bk+1][bj+1][bi+1].z + w[1]*fieldVec[bk+1][bj+1][bi+2].z
                      + w[2]*fieldVec[bk+1][bj+2][bi+1].z + w[3]*fieldVec[bk+1][bj+2][bi+2].z
                      + w[4]*fieldVec[bk+2][bj+1][bi+1].z + w[5]*fieldVec[bk+2][bj+1][bi+2].z
                      + w[6]*fieldVec[bk+2][bj+2][bi+1].z + w[7]*fieldVec[bk+2][bj+2][bi+2].z;
            } else {
                ierr = TrilinearInterpolation_Vector(fieldName, fieldVec,
                    bi + 1, bj + 1, bk + 1, a1_new, a2_new, a3_new, &vec); CHKERRQ(ierr);
            }
            ((PetscReal*)swarmOut)[3*p + 0] = vec.x;
            ((PetscReal*)swarmOut)[3*p + 1] = vec.y;
            ((PetscReal*)swarmOut)[3*p + 2] = vec.z;
        }
    }
 
    /* Restore arrays and swarm fields */
    ierr = DMDAVecRestoreArrayRead(dm_field, fieldLocal_cellCentered, &fieldPtr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCent, (void *)&cent); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished direct cell-center interpolation for field '%s'.\n", fieldName);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateEulerFieldFromCornerToSwarm"
 
/**
 * @brief Corner-averaged interpolation path (legacy).
 * @details Stages cell-centered data to corners via unweighted averaging,
 *          then trilinear interpolation from corners to particles.
 *          Local to this translation unit.
 */
static PetscErrorCode InterpolateEulerFieldFromCornerToSwarm(
    UserCtx    *user,
    Vec         fieldLocal_cellCentered,
    const char *fieldName,
    const char *swarmOutFieldName)
{
    PetscErrorCode ierr;
    DM             fda = user->fda;
    DM             swarm = user->swarm;
    PetscInt       bs;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
 
    // Generic pointers to the raw data arrays
    void *cellCenterPtr_read;
    void *cornerPtr_read_with_ghosts;
 
    // Swarm-related pointers
    PetscInt  *cellIDs = NULL;
    PetscReal *weights = NULL;
    PetscInt64  *pids = NULL; // For logging particle IDs in warnings
    void      *swarmOut = NULL;
    PetscReal *pos = NULL;
    PetscInt  *status = NULL;
    PetscInt   nLocal;
 
    PetscFunctionBegin;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(fda, &info); CHKERRQ(ierr);
    ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
    if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting for field '%s'... \n", fieldName);
    // Select the appropriate DM for intermediate corner data based on block size
    DM dm_corner = (bs == 3) ? user->fda : user->da;
 
    /* STAGE 1: Center-to-Corner Calculation with Communication */
 
    // (A) Create or reuse the corner vectors.
    if(user->CellFieldAtCorner){
        PetscInt cached_bs;
        ierr = VecGetBlockSize(user->CellFieldAtCorner, &cached_bs); CHKERRQ(ierr);
        if(cached_bs != bs){
                LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d] Block size mismatch for existing corner vector and current field %s(Cached: %d, New: %d). Replacing.\n", rank, fieldName, cached_bs, bs);
            ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr);
            ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr);
            user->CellFieldAtCorner = NULL;
            user->lCellFieldAtCorner = NULL;    
        }
    }
 
    if(user->CellFieldAtCorner == NULL){
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
        ierr = DMCreateGlobalVector(dm_corner, &user->CellFieldAtCorner); CHKERRQ(ierr);
    }else{
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing global corner vector for field '%s'.\n", rank, fieldName);
    }
 
    ierr = VecSet(user->CellFieldAtCorner,0.0); CHKERRQ(ierr);
    Vec cornerGlobal = user->CellFieldAtCorner;
 
    if(user->lCellFieldAtCorner == NULL){
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new local corner vector for field '%s'.\n", rank, fieldName);
        ierr = DMCreateLocalVector(dm_corner, &user->lCellFieldAtCorner); CHKERRQ(ierr);
    }else{
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing local corner vector for field '%s'.\n", rank, fieldName);
    }
 
    ierr = VecSet(user->lCellFieldAtCorner,0.0); CHKERRQ(ierr);
    Vec cornerLocal = user->lCellFieldAtCorner;
 
    // (B) Get a read-only array from the input cell-centered vector
    ierr = DMDAVecGetArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
 
    PetscInt size; 
    ierr = VecGetSize(cornerGlobal, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
    
    size = 0;
 
    ierr = VecGetSize(cornerLocal, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
 
    size = 0;
    
    ierr = VecGetSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
    
    PetscInt xs,ys,zs,gxs,gys,gzs;
    
    ierr  = DMDAGetCorners(dm_corner,&xs,&ys,&zs,NULL,NULL,NULL); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetCorners for field '%s': xs=%d, ys=%d, zs=%d.\n", rank, fieldName, (PetscInt)xs, (PetscInt)ys, (PetscInt)zs);
 
    ierr  = DMDAGetGhostCorners(dm_corner,&gxs,&gys,&gzs,NULL,NULL,NULL); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetGhostCorners for field '%s': gxs=%d, gys=%d, gzs=%d.\n", rank, fieldName, (PetscInt)gxs, (PetscInt)gys, (PetscInt)gzs);
    
    // DEBUG: Inspect Ucat ghost cells before interpolation
    /*
    if (bs == 3) {
        Cmpnts ***ucat_array = (Cmpnts***)cellCenterPtr_read;
        DMDALocalInfo info_debug;
        ierr = DMDAGetLocalInfo(fda, &info_debug); CHKERRQ(ierr);
 
        // Only print on rank 0 to avoid clutter
        if (rank == 0) {
            PetscPrintf(PETSC_COMM_SELF, "\nDEBUG: Inspecting Ucat Ghost Cells...\n");
            PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n");
 
            // --- Check the MINIMUM-SIDE (-Xi) ---
            // The first physical cell is at index info.xs (local index).
            // The first ghost cell is at info.xs - 1.
            PetscInt i_first_phys = info_debug.xs;
            PetscInt i_first_ghost = info_debug.xs - 1;
            PetscInt j_mid = info_debug.ys + info_debug.ym / 2;
            PetscInt k_mid = info_debug.zs + info_debug.zm / 2;
 
            PetscPrintf(PETSC_COMM_SELF, "MIN-SIDE (-Xi):\n");
            PetscPrintf(PETSC_COMM_SELF, "  Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_first_ghost,
                        ucat_array[k_mid][j_mid][i_first_ghost].x,
                        ucat_array[k_mid][j_mid][i_first_ghost].y,
                        ucat_array[k_mid][j_mid][i_first_ghost].z);
            PetscPrintf(PETSC_COMM_SELF, "  First Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_first_phys,
                        ucat_array[k_mid][j_mid][i_first_phys].x,
                        ucat_array[k_mid][j_mid][i_first_phys].y,
                        ucat_array[k_mid][j_mid][i_first_phys].z);
 
 
            // --- Check the MAXIMUM-SIDE (+Xi) ---
            // The last physical cell is at index info.xs + info.xm - 1.
            // The first ghost cell on the max side is at info.xs + info.xm.
            PetscInt i_last_phys = info_debug.xs + info_debug.xm - 1;
            PetscInt i_last_ghost = info_debug.xs + info_debug.xm;
 
            PetscPrintf(PETSC_COMM_SELF, "MAX-SIDE (+Xi):\n");
            PetscPrintf(PETSC_COMM_SELF, "  Last Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_last_phys,
                        ucat_array[k_mid][j_mid][i_last_phys].x,
                        ucat_array[k_mid][j_mid][i_last_phys].y,
                        ucat_array[k_mid][j_mid][i_last_phys].z);
            PetscPrintf(PETSC_COMM_SELF, "  Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_last_ghost,
                        ucat_array[k_mid][j_mid][i_last_ghost].x,
                        ucat_array[k_mid][j_mid][i_last_ghost].y,
                        ucat_array[k_mid][j_mid][i_last_ghost].z);
            PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n\n");
        }
    }
    */
    // DEBUG   
 
    // (C) Perform the center-to-corner interpolation directly into the global corner vector
    void *cornerPtr_write = NULL;
    ierr = DMDAVecGetArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Starting center-to-corner interpolation for '%s'.\n", rank, fieldName);
        
    // SINGLE, CLEAN CALL SITE: The macro handles the runtime dispatch based on 'bs'.
    ierr = InterpolateFieldFromCenterToCorner(bs, cellCenterPtr_read, cornerPtr_write, user); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Finished center-to-corner interpolation for '%s'.\n", rank, fieldName);
    ierr = DMDAVecRestoreArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
    
    ierr = DMDAVecRestoreArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
 
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    //////////// DEBUG
    /*
        // Log a synchronized header message from all ranks.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerGlobal BEFORE ghost exchange...\n");
 
    // VecView prints to a PETSc viewer. PETSC_VIEWER_STDOUT_WORLD is a global, synchronized viewer.
    // We only need to call it if logging is active for this function.
    if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
        ierr = VecView(cornerGlobal, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
    
    // Log a synchronized footer message.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerGlobal.\n");
    */
    ////////// DEBUG
 
    // ierr = PetscBarrier((PetscObject)cornerGlobal); CHKERRQ(ierr);
    
    // (D) CRITICAL STEP: Communicate the newly computed corner data to fill ghost regions
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Beginning ghost exchange for corner data...\n", rank);
 
    ierr = DMGlobalToLocalBegin(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
    
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Ghost exchange for corner data complete.\n", rank);
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE, "getting array from %s on all ranks.\n", dm_corner == user->fda ? "fda" : "da");
    if (is_function_allowed(__func__) && (int)(LOG_VERBOSE) <= (int)get_log_level()) {
        // DEBUG: Inspect cornerLocal after ghost exchange
        
        ierr = LOG_FIELD_ANATOMY(user,"CornerField", "After Corner Velocity Interpolated"); CHKERRQ(ierr);
 
        // DEBUG
        //ierr = DMView(user->fda, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
    
    /////// DEBUG
    /*
        // Log a synchronized header message from all ranks.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerLocal AFTER ghost exchange...\n");
 
    // Here, we want each rank to print its own local vector.
    // PETSC_VIEWER_STDOUT_SELF is the correct tool for this. The LOG_ALLOW_SYNC
    // wrapper will ensure the output from each rank is printed sequentially.
    if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
        PetscMPIInt rank_d;
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_d); CHKERRQ(ierr);
        
        // Use a synchronized print to create a clean header for each rank's output
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n--- cornerLocal (Rank %d) ---\n", rank_d);
        PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
 
        // Print the local vector's contents for this rank
        ierr = VecView(cornerLocal, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
 
        // Use a synchronized print to create a clean footer
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--- End cornerLocal (Rank %d) ---\n", rank_d);
        PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerLocal.\n");
    */
    ///// DEBUG
    
    // STAGE 2: Particle Interpolation using Ghosted Corner Data */
 
    // (E) Get the local, GHOSTED array of corner data for the final interpolation step
    ierr = DMDAVecGetArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
 
    // (F) Retrieve swarm fields for the particle loop
    ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL,LOG_TRACE," Rank %d holds data upto & including %d,%d,%d.\n",rank,info.gxs + info.gxm,info.gys+info.gym,info.gzs+info.gzm);
 
    // (G) Loop over each local particle
    for (PetscInt p = 0; p < nLocal; p++) {
        
        Particle particle;
 
        ierr = UnpackSwarmFields(p,pids,weights,pos,cellIDs,NULL,status,NULL,NULL,NULL,&particle); CHKERRQ(ierr);
        
        LOG_ALLOW(LOCAL, LOG_VERBOSE,
          "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
          rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2],
          particle.weights.x, particle.weights.y, particle.weights.z);
        // Safety check: Ensure the entire 8-node stencil is within the valid memory
        // region (owned + ghosts) of the local array.
 
        if (particle.cell[0] < info.gxs || particle.cell[0] >= info.gxs + info.gxm - 1 ||
            particle.cell[1] < info.gys || particle.cell[1] >= info.gys + info.gym - 1 ||
            particle.cell[2] < info.gzs || particle.cell[2] >= info.gzs + info.gzm - 1)
        {
            LOG_ALLOW(LOCAL, LOG_WARNING,
                      "[Rank %d] Particle PID %lld in global cell (%d,%d,%d) is in an un-interpolatable region (requires ghosts of ghosts or is out of bounds). Zeroing field '%s'.\n",
                      rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2], fieldName);
            if (bs == 3) {
                ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
            } else {
                 ((PetscReal*)swarmOut)[p] = 0.0;
            }
            continue;
        }
 
        ierr = InterpolateEulerFieldToSwarmForParticle(
                  fieldName,
                  cornerPtr_read_with_ghosts,
                  &particle,
                  swarmOut, p, bs);
        CHKERRQ(ierr);
    }
 
    // (H) Restore all retrieved arrays and vectors */
    ierr = DMDAVecRestoreArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
    
    // (I) Restore swarm fields
    ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished for field '%s'.\n", fieldName);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateEulerFieldToSwarm"
 
/**
 * @brief Dispatches grid-to-particle interpolation to the method selected in the control file.
 * @details Routes to InterpolateEulerFieldFromCenterToSwarm (direct trilinear, second-order)
 *          or InterpolateEulerFieldFromCornerToSwarm (corner-averaged, legacy) based on
 *          user->simCtx->interpolationMethod.
 */
PetscErrorCode InterpolateEulerFieldToSwarm(
    UserCtx    *user,
    Vec         fieldLocal_cellCentered,
    const char *fieldName,
    const char *swarmOutFieldName)
{
    PetscErrorCode ierr;
    PetscFunctionBegin;
 
    if (user->simCtx->interpolationMethod == INTERP_TRILINEAR) {
        ierr = InterpolateEulerFieldFromCenterToSwarm(user, fieldLocal_cellCentered,
                                                       fieldName, swarmOutFieldName); CHKERRQ(ierr);
    } else {
        ierr = InterpolateEulerFieldFromCornerToSwarm(user, fieldLocal_cellCentered,
                                                       fieldName, swarmOutFieldName); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateAllFieldsToSwarm"
/**
 * @brief Internal helper implementation: `InterpolateAllFieldsToSwarm()`.
 * @details Local to this translation unit.
 */
PetscErrorCode  InterpolateAllFieldsToSwarm(UserCtx *user)
{
  PetscErrorCode ierr;
  PetscMPIInt rank;
  PetscFunctionBegin;
 
  PROFILE_FUNCTION_BEGIN;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
 
  /* 
     1) Interpolate the 'velocity' field (user->Ucat) into DMSwarm's "swarmVelocity".
        - The function InterpolateOneFieldOverSwarm() loops over all local particles,
          retrieves iCell, jCell, kCell, (a1,a2,a3) from the DMSwarm,
          and calls the appropriate trilinear interpolation routine.
 
        - "velocity" => just a label used in logging or debugging.
        - "swarmVelocity" => the DMSwarm field name where we store the result
          (assume dof=3 if user->Ucat has blockSize=3).
  */
 
  LOG_ALLOW(LOCAL, LOG_DEBUG,
         " Interpolation of ucat to velocity begins on rank %d.\n",rank);
  // Make sure to pass the *LOCAL* Vector to the function below! 
  ierr = InterpolateEulerFieldToSwarm(user, user->lUcat, 
                                      "Ucat", 
                                      "velocity"); CHKERRQ(ierr);
  ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivity,
                                      "Diffusivity", "Diffusivity"); CHKERRQ(ierr);
  
  ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivityGradient,
                                      "DiffusivityGradient", "DiffusivityGradient"); CHKERRQ(ierr);
  /* Add fields as necessary here*/
                                      
  ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
  LOG_ALLOW(LOCAL, LOG_INFO,
        "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}
 
/////////////////////// Scatter from particles to euler fields
 
/**
 * @brief Functions for scattering particle data (scalar or vector) onto
 *        Eulerian grid fields by averaging contributions within each cell.
 *
 * This file provides a modular set of functions to perform particle-to-grid
 * projection, specifically calculating cell-averaged quantities from particle properties.
 * It assumes a PETSc environment using DMDA for the grids and DMSwarm for particles.
 *
 * Key Features:
 * - Handles both scalar (DOF=1) and vector (DOF=3) particle fields.
 * - Uses a pre-calculated particle count vector (`ParticleCount`) for normalization.
 * - Implicitly determines the target Eulerian DM (typically `user->da` for scalars,
 *   `user->fda` for vectors) based on standard field names ("P", "Nvert", "Ucat", "Ucont","Psi").
 * - Modifies existing Eulerian field vectors (e.g., `user->P`, `user->Ucat`,user->Psi`) in place.
 * - Provides a high-level wrapper function (`ScatterAllParticleFieldsToEulerFields`)
 *   to easily scatter a standard set of fields.
 * - Uses only the base `SETERRQ` macro for error reporting to maximize compiler compatibility.
 *
 * Dependencies:
 * - PETSc library (DMDA, DMSwarm, Vec, IS, PetscLog, etc.)
 * - A `UserCtx` struct (defined elsewhere, e.g., "userctx.h") containing pointers
 *   to relevant DMs (`da`, `fda`), Vecs (`ParticleCount`, `P`, `Nvert`, `Ucat`, etc.),
 *   and the `DMSwarm` object (`swarm`).
 * - A custom DMSwarm field named `"DMSwarm_CellID"` (blockSize=3, type=PETSC_INT)
 *   must be registered and populated with the local cell indices for each particle.
 * - Logging infrastructure (`LOG_ALLOW`, etc.) assumed to be defined elsewhere.
 *
 * @defgroup scatter_module Particle-to-Grid Scattering
 * @{
 */
 
//-----------------------------------------------------------------------------
// Internal Helper Modules (Lower-level building blocks)
//-----------------------------------------------------------------------------
/**
 * @defgroup scatter_module_internal Internal Scattering Helpers
 * @ingroup scatter_module
 * @brief Lower-level functions used by the main scattering routines.
 * @{
 */
 
#undef __FUNCT__
#define __FUNCT__ "GetScatterTargetInfo"
/* Implementation for GetScatterTargetInfo declared in include/interpolation.h. */
PetscErrorCode GetScatterTargetInfo(UserCtx *user, const char *particleFieldName,
                                    DM *targetDM, PetscInt *expected_dof)
{
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- Input Validation ---
    // Check for NULL pointers in essential inputs
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL.");
    if (!user->da) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->da is NULL.");
    if (!user->fda) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->fda is NULL."); // Needed for vector fields
    if (!particleFieldName) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "particleFieldName is NULL.");
    if (!targetDM) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output targetDM pointer is NULL.");
    if (!expected_dof) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output expected_dof pointer is NULL.");
 
    // --- Determine Target DM and DOF based on Field Name ---
    // Compare the input field name with known scalar fields targeting 'da'
    if (strcmp(particleFieldName, "Psi") == 0 || strcmp(particleFieldName, "Nvert") == 0) {
        *expected_dof = 1;      // Scalar fields have DOF 1
        *targetDM = user->da;   // Target the primary scalar DMDA
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Field '%s' targets DM 'da' (DOF=1).\n", particleFieldName);
    }
    // Compare with known vector fields targeting 'fda'
    else if (strcmp(particleFieldName, "Ucat") == 0 || strcmp(particleFieldName, "Ucont") == 0) {
        *expected_dof = 3;      // Vector fields have DOF 3
        *targetDM = user->fda;  // Target the vector DMDA (often node-based)
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Field '%s' targets DM 'fda' (DOF=3).\n", particleFieldName);
    }
    // --- Add rules for other fields here ---
    // else if (strcmp(particleFieldName, "SomeOtherScalar") == 0) { *expected_dof = 1; *targetDM = user->da; }
    // else if (strcmp(particleFieldName, "SomeOtherVector") == 0) { *expected_dof = 3; *targetDM = user->someOtherDM; }
    else {
        // The provided field name doesn't match any known rules
        *targetDM = NULL; // Indicate failure
        *expected_dof = 0;
        // Format the error message manually
        PetscSNPrintf(msg, sizeof(msg), "Field name '%s' is not recognized for automatic DM selection.", particleFieldName);
        // Use SETERRQ with the formatted message and appropriate error code
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg); // Use WRONG argument error code
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "GetPersistentLocalVector"
/**
 * @brief Internal helper implementation: `GetPersistentLocalVector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode GetPersistentLocalVector(UserCtx *user, const char *fieldName, Vec *localVec)
{
    PetscFunctionBeginUser;
    
    if (!user || !fieldName || !localVec) 
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null input in GetPersistentLocalVector");
 
    if (strcmp(fieldName, "Psi") == 0) {
        if (!user->lPsi) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->lPsi is not initialized");
        *localVec = user->lPsi;
    }
    else if (strcmp(fieldName, "Ucat") == 0) {
        if (!user->lUcat) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->lUcat is not initialized");
        *localVec = user->lUcat;
    }
    else if (strcmp(fieldName, "Ucont") == 0) {
        if (!user->lUcont) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->lUcont is not initialized");
        *localVec = user->lUcont;
    }
    else if (strcmp(fieldName, "Nvert") == 0) {
        if (!user->lNvert) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->lNvert is not initialized");
        *localVec = user->lNvert;
    }
    else {
        char msg[ERROR_MSG_BUFFER_SIZE];
        PetscSNPrintf(msg, sizeof(msg), "Field '%s' does not have a mapped persistent local vector.", fieldName);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg);
    }
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "AccumulateParticleField"
/*
 * Internal implementation detail for public API `AccumulateParticleField`.
 * Documentation is maintained in include/interpolation.h to keep a single source
 * of truth for Doxygen output.
 */
PetscErrorCode AccumulateParticleField(DM swarm, const char *particleFieldName,
                                       DM gridSumDM, Vec localAccumulatorVec)
{
    PetscErrorCode    ierr;
    PetscInt          dof;
    PetscInt          nlocal, p;
    const PetscReal   *particle_arr = NULL;
    const PetscInt    *cell_id_arr = NULL;
    
    // DMDA Accessors
    PetscScalar       ***arr_1d = NULL;      // For scalar fields
    PetscScalar       ****arr_3d = NULL;     // For vector fields
    
    // Ghosted dimension variables
    PetscInt          gxs, gys, gzs, gxm, gym, gzm;
    PetscMPIInt       rank;
    char              msg[ERROR_MSG_BUFFER_SIZE];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- 1. Validation & Setup ---
    if (!swarm || !gridSumDM || !localAccumulatorVec) 
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null input in AccumulateParticleField.");
 
    // Get DMDA information
    ierr = DMDAGetInfo(gridSumDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    // Get Ghosted Corners (Global indices of the ghosted patch start, and its dimensions)
    ierr = DMDAGetGhostCorners(gridSumDM, &gxs, &gys, &gzs, &gxm, &gym, &gzm); CHKERRQ(ierr);
 
    // --- 2. Verify Vector Type (Global vs Local) ---
    {
        PetscInt vecSize;
        PetscInt expectedLocalSize = gxm * gym * gzm * dof;
        ierr = VecGetSize(localAccumulatorVec, &vecSize); CHKERRQ(ierr);
        
        if (vecSize != expectedLocalSize) {
            PetscSNPrintf(msg, sizeof(msg), 
                "Vector dimension mismatch! Expected Ghosted Local Vector size %d (gxm*gym*gzm*dof), got %d. "
                "Did you pass a Global Vector instead of a Local Vector?", 
                expectedLocalSize, vecSize);
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg);
        }
    }
 
    // --- 3. Acquire Particle Data ---
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    // These calls will fail nicely if the field doesn't exist
    ierr = DMSwarmGetField(swarm, particleFieldName, NULL, NULL, (void **)&particle_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_id_arr); CHKERRQ(ierr);
 
    // --- 4. Acquire Grid Accessors ---
    // DMDAVecGetArray* handles the mapping from Global (i,j,k) to the underlying Local Array index
    if (dof == 1) {
        ierr = DMDAVecGetArray(gridSumDM, localAccumulatorVec, &arr_1d); CHKERRQ(ierr);
    } else if (dof == 3) {
        ierr = DMDAVecGetArrayDOF(gridSumDM, localAccumulatorVec, &arr_3d); CHKERRQ(ierr);
    } else {
         PetscSNPrintf(msg, sizeof(msg), "Unsupported DOF=%d. AccumulateParticleField supports DOF 1 or 3.", dof);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "%s", msg);
    }
 
    // --- 5. Accumulate Loop ---
    // Iterate through all local particles
    for (p = 0; p < nlocal; ++p) {
        // Retrieve Geometric Grid Index (0-based)
        PetscInt i_geom = cell_id_arr[p * 3 + 0];
        PetscInt j_geom = cell_id_arr[p * 3 + 1];
        PetscInt k_geom = cell_id_arr[p * 3 + 2];
 
        // Apply Shift (+1) to match Memory Layout (Index 0 is boundary/ghost)
        PetscInt i = i_geom + 1;
        PetscInt j = j_geom + 1;
        PetscInt k = k_geom + 1;
 
        // Bounds Check: Ensure (i,j,k) falls within the Local Ghosted Patch.
        // This allows writing to ghost slots which will later be reduced to the owner rank.
        if (i >= gxs && i < gxs + gxm && 
            j >= gys && j < gys + gym && 
            k >= gzs && k < gzs + gzm)
        {
            if (dof == 1) {
                arr_1d[k][j][i] += particle_arr[p];
            } else {
                // For DOF=3, unroll the loop for slight optimization
                arr_3d[k][j][i][0] += particle_arr[p * 3 + 0];
                arr_3d[k][j][i][1] += particle_arr[p * 3 + 1];
                arr_3d[k][j][i][2] += particle_arr[p * 3 + 2];
            }
        } 
        // Note: Particles outside the ghost layer are skipped. This is expected behavior
        // if particles have not yet been migrated or localized correctly.
    }
 
    // --- 6. Restore Arrays and Fields ---
    if (dof == 1) {
        ierr = DMDAVecRestoreArray(gridSumDM, localAccumulatorVec, &arr_1d); CHKERRQ(ierr);
    } else {
        ierr = DMDAVecRestoreArrayDOF(gridSumDM, localAccumulatorVec, &arr_3d); CHKERRQ(ierr);
    }
    
    ierr = DMSwarmRestoreField(swarm, particleFieldName, NULL, NULL, (void **)&particle_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_id_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "NormalizeGridVectorByCount"
 
/* Implementation for NormalizeGridVectorByCount declared in include/interpolation.h. */
PetscErrorCode NormalizeGridVectorByCount(DM countDM, Vec countVec,
                                          DM dataDM, Vec sumVec, Vec avgVec)
{
    PetscErrorCode ierr;
    PetscInt       data_dof;
    PetscInt       count_dof;
    PetscMPIInt    rank;
    char           msg[ERROR_MSG_BUFFER_SIZE];
 
    // Pointers for DMDA array accessors - declare specific types
    PetscScalar    ***count_arr_3d = NULL;      // For DOF=1 count vector (3D DMDA)
    PetscScalar    ***sum_arr_scalar = NULL;    // For DOF=1 sum vector (3D DMDA)
    PetscScalar    ***avg_arr_scalar = NULL;    // For DOF=1 avg vector (3D DMDA)
    PetscScalar    ****sum_arr_vector = NULL;   // For DOF=3 sum vector (3D DMDA + DOF)
    PetscScalar    ****avg_arr_vector = NULL;   // For DOF=3 avg vector (3D DMDA + DOF)
 
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- Validation ---
    ierr = DMDAGetInfo(countDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &count_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    ierr = DMDAGetInfo(dataDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &data_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    if (count_dof != 1) { PetscSNPrintf(msg, sizeof(msg), "countDM must have DOF=1, got %d.", count_dof); SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg); }
    if (data_dof != 1 && data_dof != 3) { PetscSNPrintf(msg, sizeof(msg), "dataDM DOF must be 1 or 3, got %d.", data_dof); SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg); }
 
    // --- Get Array Access using appropriate DMDA accessors ---
    ierr = DMDAVecGetArrayRead(countDM, countVec, &count_arr_3d); CHKERRQ(ierr);
 
    if (data_dof == 1) {
        ierr = DMDAVecGetArrayRead(dataDM, sumVec, &sum_arr_scalar); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dataDM, avgVec, &avg_arr_scalar); CHKERRQ(ierr);
    } else { // data_dof == 3
        ierr = DMDAVecGetArrayDOFRead(dataDM, sumVec, &sum_arr_vector); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayDOF(dataDM, avgVec, &avg_arr_vector); CHKERRQ(ierr);
    }
 
    // Get the corners (global start indices) and dimensions of the *local owned* region
    PetscInt xs, ys, zs, xm, ym, zm;
    ierr = DMDAGetCorners(countDM, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ(ierr);
 
    // --- Normalize Over Owned Cells ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "(Rank %d): Normalizing DOF=%d data over owned range [%d:%d, %d:%d, %d:%d].\n",
              rank, data_dof, xs, xs+xm, ys, ys+ym, zs, zs+zm);
 
    // Loop using GLOBAL indices (i, j, k) over the range owned by this process
    for (PetscInt k = zs; k < zs + zm; ++k) {
        for (PetscInt j = ys; j < ys + ym; ++j) {
            for (PetscInt i = xs; i < xs + xm; ++i) {
 
                // Access the count using standard 3D indexing
                PetscScalar count = count_arr_3d[k][j][i];
 
                if (PetscRealPart(count) > 0.5) { // Use tolerance for float comparison
                    if (data_dof == 1) {
                        // Access scalar sum/avg using standard 3D indexing
                        avg_arr_scalar[k][j][i] = sum_arr_scalar[k][j][i] / count;
                    } else { // data_dof == 3
                        // Access vector components using DOF indexing on the last dimension
                        for (PetscInt c = 0; c < data_dof; ++c) {
                            avg_arr_vector[k][j][i][c] = sum_arr_vector[k][j][i][c] / count;
                        }
                    }
                } else { // count is zero or negative
                    // Set average to zero
                    if (data_dof == 1) {
                        avg_arr_scalar[k][j][i] = 0.0;
                    } else { // data_dof == 3
                        for (PetscInt c = 0; c < data_dof; ++c) {
                            avg_arr_vector[k][j][i][c] = 0.0;
                        }
                    }
                } // end if count > 0.5
            } // end i loop
        } // end j loop
    } // end k loop
 
    // --- Restore Arrays using appropriate functions ---
    ierr = DMDAVecRestoreArrayRead(countDM, countVec, &count_arr_3d); CHKERRQ(ierr);
    if (data_dof == 1) {
        ierr = DMDAVecRestoreArrayRead(dataDM, sumVec, &sum_arr_scalar); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dataDM, avgVec, &avg_arr_scalar); CHKERRQ(ierr);
    } else { // data_dof == 3
        ierr = DMDAVecRestoreArrayDOFRead(dataDM, sumVec, &sum_arr_vector); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayDOF(dataDM, avgVec, &avg_arr_vector); CHKERRQ(ierr);
    }
 
    // --- Assemble Final Average Vector ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling final average vector (DOF=%d).\n", data_dof);
    ierr = VecAssemblyBegin(avgVec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(avgVec); CHKERRQ(ierr);
 
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
/** @} */ // End of scatter_module_internal group
 
//-----------------------------------------------------------------------------
// User-Facing API
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
// MODULE 4: Internal Scatter Orchestration Helper - No PetscErrorClear
//-----------------------------------------------------------------------------
 
#undef __FUNCT__
#define __FUNCT__   "ScatterParticleFieldToEulerField_Internal"
/**
 * @brief Internal helper implementation: `ScatterParticleFieldToEulerField_Internal()`.
 * @details Local to this translation unit.
 */
static PetscErrorCode ScatterParticleFieldToEulerField_Internal(UserCtx *user,
                                                         const char *particleFieldName,
                                                         DM targetDM,
                                                         PetscInt expected_dof,
                                                         Vec eulerFieldAverageVec)
{
    PetscErrorCode ierr;
    PetscInt       target_dof = 0;
    Vec            globalsumVec = NULL;
    Vec            localsumVec = NULL;
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (!user || !user->swarm || !user->ParticleCount || !particleFieldName || !targetDM || !eulerFieldAverageVec)
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "NULL input provided to ScatterParticleFieldToEulerField_Internal.");
 
    ierr = DMDAGetInfo(targetDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &target_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    if (target_dof != expected_dof) {
        PetscSNPrintf(msg, sizeof(msg),
                      "Field '%s' expects DOF %d but targetDM reports DOF %d.",
                      particleFieldName, expected_dof, target_dof);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "%s", msg);
    }
    
    // --- Check if Particle Field Exists ---
    // Attempt a GetField call; if it fails, the field doesn't exist.
    // We let CHKERRQ handle the error directly if the field doesn't exist OR
    // we catch it specifically to provide a more tailored message.
 
    /*
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"Field %s being accessed to check existence \n",particleFieldName);
    ierr = DMSwarmGetField(user->swarm, particleFieldName, NULL, NULL, NULL);
    if (ierr) { // If GetField returns an error
        PetscSNPrintf(msg, sizeof(msg), "Particle field '%s' not found in DMSwarm for scattering.", particleFieldName);
        // Directly set the error, overwriting the one from GetField
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, msg);
    }
    ierr = DMSwarmRestoreField(user->swarm, particleFieldName, NULL, NULL, NULL);
    */
    
    // --- Setup Temporary Sum Vector ---
    ierr = VecDuplicate(eulerFieldAverageVec, &globalsumVec); CHKERRQ(ierr);
    ierr = VecSet(globalsumVec, 0.0); CHKERRQ(ierr);
    ierr = PetscSNPrintf(msg, sizeof(msg), "TempSum_%s", particleFieldName); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)globalsumVec, msg); CHKERRQ(ierr);
 
    // create local vector for accumulation
    ierr = DMGetLocalVector(targetDM, &localsumVec); CHKERRQ(ierr);
    ierr = VecSet(localsumVec, 0.0); CHKERRQ(ierr);  // Must be zeroed before accumulation
    ierr = PetscSNPrintf(msg, sizeof(msg), "LocalTempSum_%s", particleFieldName); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)localsumVec, msg); CHKERRQ(ierr);
 
    // --- Accumulate ---
    // This will call DMSwarmGetField again. If it failed above, it will likely fail here too,
    // unless the error was cleared somehow between the check and here (unlikely).
    // If the check above was skipped (Option 1), this is where the error for non-existent
    // field will be caught by CHKERRQ.
    ierr = AccumulateParticleField(user->swarm, particleFieldName, targetDM, localsumVec); CHKERRQ(ierr);
 
    // --- Local to Global Sum ---
    ierr = DMLocalToGlobalBegin(targetDM, localsumVec, ADD_VALUES, globalsumVec); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(targetDM, localsumVec, ADD_VALUES, globalsumVec); CHKERRQ(ierr);
    // Return local vector to DM
    ierr = DMRestoreLocalVector(targetDM, &localsumVec); CHKERRQ(ierr);
 
    // Calculate the number of particles per cell.
    ierr = CalculateParticleCountPerCell(user); CHKERRQ(ierr);
    // --- Normalize ---
    ierr = NormalizeGridVectorByCount(user->da, user->ParticleCount, targetDM, globalsumVec, eulerFieldAverageVec); CHKERRQ(ierr);
 
    // --- Cleanup ---
    ierr = VecDestroy(&globalsumVec); CHKERRQ(ierr);
 
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ScatterParticleFieldToEulerField"
 
/* Implementation for ScatterParticleFieldToEulerField declared in include/interpolation.h. */
PetscErrorCode ScatterParticleFieldToEulerField(UserCtx *user,
                                                const char *particleFieldName,
                                                Vec eulerFieldAverageVec)
{
    PetscErrorCode ierr;
    DM             targetDM = NULL;       // Will point to user->da or user->fda
    PetscInt       expected_dof = 0;      // Will be 1 or 3
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- Essential Input Validation ---
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL.");
    if (!user->swarm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->swarm is NULL.");
    if (!user->ParticleCount) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->ParticleCount is NULL.");
    if (!eulerFieldAverageVec) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output eulerFieldAverageVec is NULL.");
    // particleFieldName validity checked within GetScatterTargetInfo
 
    // --- Determine Target DM & DOF using the helper function ---
    ierr = GetScatterTargetInfo(user, particleFieldName, &targetDM, &expected_dof); CHKERRQ(ierr);
    // If GetScatterTargetInfo failed (e.g., unknown name), CHKERRQ would have returned.
 
    // --- Validate the provided Target Vec's Compatibility ---
    DM vec_dm;
    PetscInt vec_dof;
    // Check that the provided average vector has a DM associated with it
    ierr = VecGetDM(eulerFieldAverageVec, &vec_dm); CHKERRQ(ierr);
    if (!vec_dm) {
         PetscSNPrintf(msg, sizeof(msg), "Provided eulerFieldAverageVec for field '%s' does not have an associated DM.", particleFieldName);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "%s", msg);
    }
    // Get the block size (DOF) of the provided vector
    ierr = VecGetBlockSize(eulerFieldAverageVec, &vec_dof); CHKERRQ(ierr);
    // Compare the vector's associated DM with the one determined by the field name
    if (vec_dm != targetDM) {
        const char *target_dm_name = "targetDM", *vec_dm_name = "vec_dm";
        // Get actual names if possible for a more informative error message
        PetscObjectGetName((PetscObject)targetDM, &target_dm_name);
        PetscObjectGetName((PetscObject)vec_dm, &vec_dm_name);
        PetscSNPrintf(msg, sizeof(msg), "Provided eulerFieldAverageVec associated with DM '%s', but field '%s' requires scatter to DM '%s'.", vec_dm_name, particleFieldName, target_dm_name);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "%s", msg);
    }
    // Compare the vector's DOF with the one expected for the field name
    if (vec_dof != expected_dof) {
         PetscSNPrintf(msg, sizeof(msg), "Field '%s' requires DOF %d, but provided eulerFieldAverageVec has DOF %d.", particleFieldName, expected_dof, vec_dof);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "%s", msg);
    }
 
    // --- Perform Scatter using Internal Helper ---
    // Log intent before calling the core logic
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering field '%s' (DOF=%d).\n", particleFieldName, expected_dof);
    ierr = ScatterParticleFieldToEulerField_Internal(user, // Pass user context
                                                     particleFieldName, // Name of particle field
                                                     targetDM, // Determined target DM (da or fda)
                                                     expected_dof, // Determined DOF (1 or 3)
                                                     eulerFieldAverageVec); // The output vector
    CHKERRQ(ierr); // Handle potential errors from the internal function
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Successfully scattered field '%s'.\n", particleFieldName);
   
    PROFILE_FUNCTION_END;
   
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "ScatterAllParticleFieldsToEulerFields"
/* Implementation for ScatterAllParticleFieldsToEulerFields declared in include/interpolation.h. */
PetscErrorCode ScatterAllParticleFieldsToEulerFields(UserCtx *user)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting scattering of specified particle fields to Eulerian grids.\n");
 
    // --- Pre-computation Check: Ensure Particle Counts are Ready ---
    if (!user->ParticleCount) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "UserCtx->ParticleCount is NULL. Compute counts before calling ScatterAllParticleFieldsToEulerFields.");
    }
 
    // --- Scatter Particle Field "Psi" -> Eulerian Field user->Psi (on da) ---
    // Check if the target Eulerian vector 'user->Psi' exists.
    if (user->Psi) {
        
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Psi' to user->Psi.\n");
        // Zero the target vector before accumulating the new average for this step/call.
 
    // Debug Verification ------------------------------------------------
    Vec swarm_Psi;
    PetscReal Avg_Psi,Avg_swarm_Psi;
         
    ierr = VecMean(user->Psi,&Avg_Psi);
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Scalar(Psi) before scatter: %.4f.\n",Avg_Psi);
         
    ierr = DMSwarmCreateGlobalVectorFromField(user->swarm,"Psi",&swarm_Psi);
    ierr = VecMean(swarm_Psi,&Avg_swarm_Psi);
 
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Particle Scalar(Psi): %.4f.\n",Avg_swarm_Psi);
 
    ierr = DMSwarmDestroyGlobalVectorFromField(user->swarm,"Psi",&swarm_Psi);
    // Debug----------------------------------------------------------------
      
    //ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
        // Call the unified scatter function. It will handle DM determination and validation.
        // It will also error out if the *particle* field "Psi" doesn't exist in the swarm.
        ierr = ScatterParticleFieldToEulerField(user, "Psi", user->Psi); CHKERRQ(ierr);
    ierr = VecMean(user->Psi,&Avg_Psi);
    
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Scalar(Psi) after  scatter: %.4f.\n",Avg_Psi);
    } else {
        // Only log a warning if the target Eulerian field is missing in the context.
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Psi': UserCtx->Psi is NULL.\n");
    }
 
    // --- (Commented Out) Scatter Other Fields ---
    // To enable scattering for other fields, uncomment the relevant block
    // AND ensure:
    // 1. The corresponding particle field (e.g., "Nvert", "Ucat") exists in user->swarm.
    // 2. The corresponding target Eulerian Vec (e.g., user->Nvert, user->Ucat) exists in user->ctx.
    // 3. The target Eulerian Vec is associated with the correct DM (da for Nvert, fda for Ucat).
 
    /*
    // Scatter Particle Field "Nvert" -> Eulerian Field user->Nvert (on da)
    if (user->Nvert) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Nvert' to user->Nvert.\n");
        ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Nvert", user->Nvert); CHKERRQ(ierr);
    } else {
         LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Nvert': UserCtx->Nvert is NULL.\n");
    }
    */
 
    /*
    // Scatter Particle Field "Ucat" -> Eulerian Field user->Ucat (on fda)
    if (user->Ucat) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Ucat' to user->Ucat.\n");
        ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Ucat", user->Ucat); CHKERRQ(ierr);
    } else {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Ucat': UserCtx->Ucat is NULL.\n");
    }
    */
 
    /*
     // Scatter Particle Field "Ucont" -> Eulerian Field user->Ucont (on fda)
    if (user->Ucont) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Ucont' to user->Ucont.\n");
        ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Ucont", user->Ucont); CHKERRQ(ierr);
    } else {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Ucont': UserCtx->Ucont is NULL.\n");
    }
    */
 
    // Add more fields as needed following the pattern above...
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished scattering specified particle fields.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}
 
/** @} */ // End of scatter_module group
 
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateCornerToFaceCenter_Scalar"
 
/**
 * @brief Internal helper implementation: `InterpolateCornerToFaceCenter_Scalar()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateCornerToFaceCenter_Scalar(
    PetscReal ***corner_arr,
    PetscReal ***faceX_arr,
    PetscReal ***faceY_arr,
    PetscReal ***faceZ_arr,
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo info;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Determine owned-cell ranges based on corner-node ownership
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    // Global exclusive end indices for cells
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs; // 0..xm
                // Average the four corners of the Y-Z face at X = i
                PetscReal sum = corner_arr[k  ][j  ][i]
                              + corner_arr[k+1][j  ][i]
                              + corner_arr[k  ][j+1][i]
                              + corner_arr[k+1][j+1][i];
                faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys; // 0..ym
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Z face at Y = j
                PetscReal sum = corner_arr[k  ][j][i  ]
                              + corner_arr[k+1][j][i  ]
                              + corner_arr[k  ][j][i+1]
                              + corner_arr[k+1][j][i+1];
                faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Y face at Z = k
                PetscReal sum = corner_arr[k][j  ][i  ]
                              + corner_arr[k][j  ][i+1]
                              + corner_arr[k][j+1][i  ]
                              + corner_arr[k][j+1][i+1];
                faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InterpolateCornerToFaceCenter_Vector"
 
/**
 * @brief Internal helper implementation: `InterpolateCornerToFaceCenter_Vector()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InterpolateCornerToFaceCenter_Vector(
    Cmpnts ***corner_arr,
    Cmpnts ***faceX_arr,
    Cmpnts ***faceY_arr,
    Cmpnts ***faceZ_arr,
    UserCtx *user)
{
    PetscErrorCode ierr;
    DMDALocalInfo info;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
    
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // X-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j  ][i].x + corner_arr[k+1][j  ][i].x
                      + corner_arr[k  ][j+1][i].x + corner_arr[k+1][j+1][i].x;
                sum.y = corner_arr[k  ][j  ][i].y + corner_arr[k+1][j  ][i].y
                      + corner_arr[k  ][j+1][i].y + corner_arr[k+1][j+1][i].y;
                sum.z = corner_arr[k  ][j  ][i].z + corner_arr[k+1][j  ][i].z
                      + corner_arr[k  ][j+1][i].z + corner_arr[k+1][j+1][i].z;
                faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d x-face Interpolation complete.\n", rank);
    
    // Y-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j][i  ].x + corner_arr[k+1][j][i  ].x
                      + corner_arr[k  ][j][i+1].x + corner_arr[k+1][j][i+1].x;
                sum.y = corner_arr[k  ][j][i  ].y + corner_arr[k+1][j][i  ].y
                      + corner_arr[k  ][j][i+1].y + corner_arr[k+1][j][i+1].y;
                sum.z = corner_arr[k  ][j][i  ].z + corner_arr[k+1][j][i  ].z
                      + corner_arr[k  ][j][i+1].z + corner_arr[k+1][j][i+1].z;
                faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d y-face Interpolation complete.\n", rank);
    
    // Z-faces
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k][j  ][i  ].x + corner_arr[k][j  ][i+1].x
                      + corner_arr[k][j+1][i  ].x + corner_arr[k][j+1][i+1].x;
                sum.y = corner_arr[k][j  ][i  ].y + corner_arr[k][j  ][i+1].y
                      + corner_arr[k][j+1][i  ].y + corner_arr[k][j+1][i+1].y;
                sum.z = corner_arr[k][j  ][i  ].z + corner_arr[k][j  ][i+1].z
                      + corner_arr[k][j+1][i  ].z + corner_arr[k][j+1][i+1].z;
                faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d z-face Interpolation complete.\n", rank);
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}