poisson.c

Go to the documentation of this file.

// In src/poisson.c (or wherever PoissonSolver_MG is)
#include "poisson.h" // The new header for this file
#include "logging.h"
 
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user) \
  if ((user->isc)) { \
    ic = i; \
    ia = 0; \
  } \
  else { \
    ic = (i+1) / 2; \
    ia = (i - 2 * (ic)) == 0 ? 1 : -1; \
    if (i==1 || i==mx-2) ia = 0; \
  }\
  if ((user->jsc)) { \
    jc = j; \
    ja = 0; \
  } \
  else { \
    jc = (j+1) / 2; \
    ja = (j - 2 * (jc)) == 0 ? 1 : -1; \
    if (j==1 || j==my-2) ja = 0; \
  } \
  if ((user->ksc)) { \
    kc = k; \
    ka = 0; \
  } \
  else { \
    kc = (k+1) / 2; \
    ka = (k - 2 * (kc)) == 0 ? 1 : -1; \
    if (k==1 || k==mz-2) ka = 0; \
  } \
  if (ka==-1 && nvert_c[kc-1][jc][ic] > 0.1) ka = 0; \
  else if (ka==1 && nvert_c[kc+1][jc][ic] > 0.1) ka = 0; \
  if (ja==-1 && nvert_c[kc][jc-1][ic] > 0.1) ja = 0; \
  else if (ja==1 && nvert_c[kc][jc+1][ic] > 0.1) ja = 0; \
  if (ia==-1 && nvert_c[kc][jc][ic-1] > 0.1) ia = 0; \
  else if (ia==1 && nvert_c[kc][jc][ic+1] > 0.1) ia = 0;
 
/**
 * @brief Convert local logical indices to the corresponding flattened global cell identifier.
 */
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
 
{
  PetscInt nidx;
  DMDALocalInfo info = user->info;
 
  PetscInt  mx = info.mx, my = info.my;
  
  AO ao;
  DMDAGetAO(user->da, &ao);
  nidx=i+j*mx+k*mx*my;
  
  AOApplicationToPetsc(ao,1,&nidx);
  
  return (nidx);
}
 
#undef __FUNCT__
#define __FUNCT__ "GridRestriction"
 
/**
 * @brief Restrict a fine-grid scalar value onto its associated coarse-grid location.
 */
static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k,
                   PetscInt *ih, PetscInt *jh, PetscInt *kh,
                   UserCtx *user)
{
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
  if ((user->isc)) {
    *ih = i;
  }
  else {
    *ih = 2 * i;
  }
 
  if ((user->jsc)) {
    *jh = j;
  }
  else {
    *jh = 2 * j;
  }
 
  if ((user->ksc)) {
    *kh = k;
  }
  else {
    *kh = 2 * k;
  }
 
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}
 
#include "solvers.h" // Or your main project header
 
#undef __FUNCT__
#define __FUNCT__ "CorrectChannelFluxProfile"
/**
 * @brief Internal helper implementation: `CorrectChannelFluxProfile()`.
 * @details Local to this translation unit.
 */
PetscErrorCode CorrectChannelFluxProfile(UserCtx *user)
{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
 
    // --- Step 1: Discover if and where a driven flow is active ---
    char drivenDirection = ' ';
    for (int i = 0; i < 6; i++) {
        BCHandlerType handler_type = user->boundary_faces[i].handler_type;
        if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
            handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
        {
            switch (user->boundary_faces[i].face_id) {
                case BC_FACE_NEG_X: case BC_FACE_POS_X: drivenDirection = 'X'; break;
                case BC_FACE_NEG_Y: case BC_FACE_POS_Y: drivenDirection = 'Y'; break;
                case BC_FACE_NEG_Z: case BC_FACE_POS_Z: drivenDirection = 'Z'; break;
            }
            break;
        }
    }
 
    // --- Step 2: Early exit if no driven flow is configured ---
    if (drivenDirection == ' ') {
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Starting channel flux profile correction in '%c' direction...\n",
              simCtx->rank, user->_this, drivenDirection);
 
    // --- Step 3: Setup and Get PETSc Array Pointers ---
    DMDALocalInfo info = user->info;
    PetscInt i, j, k;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
    PetscInt lys = (info.ys == 0) ? 1 : info.ys;
    PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
    PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
    PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
    PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
 
    Cmpnts ***ucont, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // --- Step 4: Allocate Memory for Profile Arrays based on direction ---
    PetscInt n_planes = 0;
    switch (drivenDirection) {
        case 'X': n_planes = mx - 1; break;
        case 'Y': n_planes = my - 1; break;
        case 'Z': n_planes = mz - 1; break;
    }
 
    PetscReal *localFluxProfile, *globalFluxProfile, *correctionProfile;
    ierr = PetscMalloc1(n_planes, &localFluxProfile); CHKERRQ(ierr);
    ierr = PetscMalloc1(n_planes, &globalFluxProfile); CHKERRQ(ierr);
    ierr = PetscMalloc1(n_planes, &correctionProfile); CHKERRQ(ierr);
    ierr = PetscMemzero(localFluxProfile, n_planes * sizeof(PetscReal)); CHKERRQ(ierr);
 
    // --- Step 5: Calculate Total Cross-Sectional Area and Measure Flux Profile ---
    PetscReal localArea = 0.0, globalArea = 0.0;
    
    switch (drivenDirection) {
        case 'X':
            if (info.xs == 0) { // Area is calculated by rank(s) on the negative face
                i = 0;
                for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
                    if (nvert[k][j][i + 1] < 0.1)
                        localArea += sqrt(csi[k][j][i].x*csi[k][j][i].x + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].z*csi[k][j][i].z);
                }
            }
            for (i = info.xs; i < lxe; i++) {
                for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
                    if (nvert[k][j][i + 1] < 0.1) localFluxProfile[i] += ucont[k][j][i].x;
                }
            }
            break;
        case 'Y':
            if (info.ys == 0) {
                j = 0;
                for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k][j + 1][i] < 0.1)
                        localArea += sqrt(eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z);
                }
            }
            for (j = info.ys; j < lye; j++) {
                for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k][j + 1][i] < 0.1) localFluxProfile[j] += ucont[k][j][i].y;
                }
            }
            break;
        case 'Z':
            if (info.zs == 0) {
                k = 0;
                for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k + 1][j][i] < 0.1)
                        localArea += sqrt(zet[k][j][i].x*zet[k][j][i].x + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].z*zet[k][j][i].z);
                }
            }
            for (k = info.zs; k < lze; k++) {
                for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k + 1][j][i] < 0.1) localFluxProfile[k] += ucont[k][j][i].z;
                }
            }
            break;
    }
 
    ierr = MPI_Allreduce(&localArea, &globalArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
    ierr = MPI_Allreduce(localFluxProfile, globalFluxProfile, n_planes, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    // --- Step 6: Calculate Correction Profile ---
    PetscReal targetFlux = simCtx->targetVolumetricFlux;
    if (globalArea > 1.0e-12) {
        for (i = 0; i < n_planes; i++) {
            correctionProfile[i] = (targetFlux - globalFluxProfile[i]) / globalArea;
        }
    } else {
        ierr = PetscMemzero(correctionProfile, n_planes * sizeof(PetscReal)); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Channel Flux Profile Corrector Update (Dir %c):\n", drivenDirection);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Target Flux for all planes: %.6e\n", targetFlux);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Measured Flux at plane 0:   %.6e (Correction Velocity: %.6e)\n", globalFluxProfile[0], correctionProfile[0]);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Measured Flux at plane %d:  %.6e (Correction Velocity: %.6e)\n", (n_planes-1)/2, globalFluxProfile[(n_planes-1)/2], correctionProfile[(n_planes-1)/2]);
 
    /*  TURNED OFF IN LEGACY 
    // --- Step 7: Apply Correction to Velocity Profile ---
    switch (drivenDirection) {
        case 'X':
            for (i = info.xs; i < info.xs + info.xm - 1; i++) {
                if (PetscAbs(correctionProfile[i]) > 1e-12) {
                    for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
                        if (nvert[k][j][i] < 0.1) {
                            PetscReal faceArea = sqrt(csi[k][j][i].x*csi[k][j][i].x + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].z*csi[k][j][i].z);
                            ucont[k][j][i].x += correctionProfile[i] * faceArea;
                        }
                    }
                }
            }
            break;
        case 'Y':
            for (j = info.ys; j < info.ys + info.ym - 1; j++) {
                if (PetscAbs(correctionProfile[j]) > 1e-12) {
                    for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
                        if (nvert[k][j][i] < 0.1) {
                            PetscReal faceArea = sqrt(eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z);
                            ucont[k][j][i].y += correctionProfile[j] * faceArea;
                        }
                    }
                }
            }
            break;
        case 'Z':
            for (k = info.zs; k < info.zs + info.zm - 1; k++) {
                if (PetscAbs(correctionProfile[k]) > 1e-12) {
                    for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
                        if (nvert[k][j][i] < 0.1) {
                            PetscReal faceArea = sqrt(zet[k][j][i].x*zet[k][j][i].x + zet[k][j][i].y*zet[k][j][i].y + zet[k][j][i].z*zet[k][j][i].z);
                            ucont[k][j][i].z += correctionProfile[k] * faceArea;
                        }
                    }
                }
            }
            break;
    }
    */
 
    // --- Step 8: Cleanup and Restore ---
    ierr = PetscFree(localFluxProfile); CHKERRQ(ierr);
    ierr = PetscFree(globalFluxProfile); CHKERRQ(ierr);
    ierr = PetscFree(correctionProfile); CHKERRQ(ierr);
 
    ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    //LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Channel flux profile correction complete.\n",
    //          simCtx->rank, user->_this);
              
    PetscFunctionReturn(0);
}
 
//Cell Center
#define CP  0 
// Face Centers
#define EP  1
#define WP  2
#define NP  3
#define SP  4
#define TP  5
#define BP  6
 
// Edge Centers
#define NE  7
#define SE  8
#define NW  9
#define SW  10
#define TN  11
#define BN  12
#define TS  13
#define BS  14
#define TE  15
#define BE  16
#define TW  17
#define BW  18
 
#undef __FUNCT__
#define __FUNCT__ "Projection"
/**
 * @brief Implementation of \ref Projection().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see Projection()
 */
PetscErrorCode Projection(UserCtx *user)
{
  PetscErrorCode ierr;
 
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;  
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering Projection step to correct velocity field.\n");
 
  //================================================================================
  // Section 1: Initialization and Data Acquisition
  //================================================================================
 
  // --- Get simulation and grid context ---
  SimCtx *simCtx = user->simCtx;
  DM da = user->da, fda = user->fda;
  DMDALocalInfo info = user->info;
 
  // --- Grid dimensions ---
  PetscInt mx = info.mx, my = info.my, mz = info.mz;
  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;
 
  // --- Loop bounds (excluding outer ghost layers) ---
  PetscInt lxs = (xs == 0) ? xs + 1 : xs;
  PetscInt lxe = (xe == mx) ? xe - 1 : xe;
  PetscInt lys = (ys == 0) ? ys + 1 : ys;
  PetscInt lye = (ye == my) ? ye - 1 : ye;
  PetscInt lzs = (zs == 0) ? zs + 1 : zs;
  PetscInt lze = (ze == mz) ? ze - 1 : ze;
 
  // --- Get direct pointer access to grid metric and field data ---
  Cmpnts ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
  PetscReal ***iaj, ***jaj, ***kaj, ***p, ***nvert;
  Cmpnts ***ucont;
  DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
  DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
  DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
  DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
  DMDAVecGetArray(da, user->lNvert, &nvert);
  DMDAVecGetArray(da, user->lPhi, &p); // Note: using lPhi, which is the pressure correction
  //DMDAVecGetArray(da,user->lP,&p);
  DMDAVecGetArray(fda, user->Ucont, &ucont);
 
  // --- Constants for clarity ---
  const PetscReal IBM_FLUID_THRESHOLD = 0.1;
  const PetscReal scale = simCtx->dt * 1.0 / COEF_TIME_ACCURACY; // simCtx->st replaced by 1.0.
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG," Starting velocity correction: Scale = %le .\n",scale);
 
  //================================================================================
  // Section 2: Correct Velocity Components
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating pressure gradients and correcting velocity components.\n");
  
  // --- Main loop over interior domain points ---
  for (PetscInt k = lzs; k < lze; k++) {
    for (PetscInt j = lys; j < lye; j++) {
      for (PetscInt i = lxs; i < lxe; i++) {
 
        // --- Correct U_contravariant (x-component of velocity) ---
        PetscInt i_end = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) ? mx - 1 : mx - 2;
        if (i < i_end) {
      
          if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j][i + 1] > IBM_FLUID_THRESHOLD)) {
            // Compute pressure derivatives (dp/d_csi, dp/d_eta, dp/d_zet) at the i-face
 
        PetscReal dpdc = p[k][j][i + 1] - p[k][j][i];
            PetscReal dpde = 0.0, dpdz = 0.0;
 
            // Boundary-aware stencil for dp/d_eta
          if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
          dpde = (p[k][j][i] + p[k][j][i+1] -
              p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
        }
          }
 
          else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) { dpde = (p[k][j][i] + p[k][j][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.5; }
          }
 
          else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
          }
 
          else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
          }
 
          else { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.25; }
 
            // Boundary-aware stencil for dp/d_zet
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
          }
 
          else  if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
          }
 
          else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
          }
 
          else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
          }
 
          else { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.25; }
 
            // Apply the correction: U_new = U_old - dt * (g11*dpdc + g12*dpde + g13*dpdz)
 
            
          
            PetscReal grad_p_x = (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y
                      + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
                                dpde * (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y
                    + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
                                dpdz * (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y
                    + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i]);
 
        PetscReal correction = grad_p_x*scale;
        //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Csi Direction: %le.\n",correction);
        ucont[k][j][i].x -= correction;
 
          }
        }
 
        // --- Correct V_contravariant (y-component of velocity) ---
        PetscInt j_end = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) ? my - 1 : my - 2;
        if (j < j_end) {
          if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j + 1][i] > IBM_FLUID_THRESHOLD)) {
            PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
            dpde = p[k][j + 1][i] - p[k][j][i];
 
            // Boundary-aware stencil for dp/d_csi
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
          } else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
          } else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
          } else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
          } else { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.25; }
 
            // Boundary-aware stencil for dp/d_zet
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
          } else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
          } else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
          } else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
          } else { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.25; }
 
            PetscReal grad_p_y = (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
                                dpde * (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
                                dpdz * (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i]);
 
        PetscReal correction = grad_p_y*scale;
        //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Eta Direction: %le.\n",correction);
            ucont[k][j][i].y -= correction;
          }
        }
        
        // --- Correct W_contravariant (z-component of velocity) ---
        PetscInt k_end = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) ? mz - 1 : mz - 2;
        if (k < k_end) {
          if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k + 1][j][i] > IBM_FLUID_THRESHOLD)) {
            PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
            dpdz = p[k + 1][j][i] - p[k][j][i];
            
            // Boundary-aware stencil for dp/d_csi
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
          } else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
          } else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
          } else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
          } else { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.25; }
 
            // Boundary-aware stencil for dp/d_eta
          if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
          } else  if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
          } else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
          } else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
          } else { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.25; }
 
            PetscReal grad_p_z = (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
                                dpde * (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
                                dpdz * (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i]);
 
        // ========================= DEBUG PRINT  =========================
        LOG_LOOP_ALLOW_EXACT(GLOBAL, LOG_DEBUG, k, 5,
                 "[k=%d, j=%d, i=%d] ---- Neighbor Pressures ----\n"
                 "  Central Z-Neighbors: p[k+1][j][i] = %g | p[k][j][i] = %g\n"
                 "  Eta-Stencil (Y-dir): p[k][j-1][i] = %g, p[k+1][j-1][i] = %g | p[k][j+1][i] = %g, p[k+1][j+1][i] = %g\n"
                 "  Csi-Stencil (X-dir): p[k][j][i-1] = %g, p[k+1][j][i-1] = %g | p[k][j][i+1] = %g, p[k+1][j][i+1] = %g\n",
                 k, j, i,
                 p[k + 1][j][i], p[k][j][i],
                 p[k][j - 1][i], p[k + 1][j - 1][i], p[k][j + 1][i], p[k + 1][j + 1][i],
                 p[k][j][i - 1], p[k + 1][j][i - 1], p[k][j][i + 1], p[k + 1][j][i + 1]);
        // ======================= END DEBUG PRINT =======================
 
        LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," dpdc: %le | dpde: %le | dpdz: %le.\n",dpdc,dpde,dpdz);      
        PetscReal correction = grad_p_z*scale;
        //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Zet Direction: %le.\n",correction);
            ucont[k][j][i].z -= correction;
          }
        }
      }
    }
  }
 
  // --- Explicit correction for periodic boundaries (if necessary) ---
  // The main loop handles the interior, but this handles the first physical layer at periodic boundaries.
  // Note: This logic is largely duplicated from the main loop and could be merged, but is preserved for fidelity.
  if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
    for (PetscInt k=lzs; k<lze; k++) {
      for (PetscInt j=lys; j<lye; j++) {
    PetscInt i=xs;
    
    PetscReal dpdc = p[k][j][i+1] - p[k][j][i];
    
    PetscReal dpde = 0.;
    PetscReal dpdz = 0.;
        
    if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
      if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
        dpde = (p[k][j  ][i] + p[k][j  ][i+1] -
            p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
      }
    }
    else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
      if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
        dpde = (p[k][j  ][i] + p[k][j  ][i+1] -
            p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
      }
    }
    else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
      if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
        dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
            p[k][j  ][i] - p[k][j  ][i+1]) * 0.5;
      }
    }
    else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
      if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
        dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
            p[k][j  ][i] - p[k][j  ][i+1]) * 0.5;
      }
    }
    else {
      dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
          p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
    }
    
    if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
      if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
        dpdz = (p[k  ][j][i] + p[k  ][j][i+1] -
            p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
      }
    }
    else  if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
      if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
        dpdz = (p[k  ][j][i] + p[k  ][j][i+1] -
            p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
      }
    }
    else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
      if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
        dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
            p[k  ][j][i] - p[k  ][j][i+1]) * 0.5;
      }
    }
    else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
      if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
        dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
            p[k  ][j][i] - p[k  ][j][i+1]) * 0.5;
      }
    }
    else {
      dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
          p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
    }
    
    
    
    if (!(nvert[k][j][i] + nvert[k][j][i+1])) {
      ucont[k][j][i].x -=
        (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
             icsi[k][j][i].y * icsi[k][j][i].y +
             icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
         dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
             ieta[k][j][i].y * icsi[k][j][i].y +
             ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
         dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
             izet[k][j][i].y * icsi[k][j][i].y +
             izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i])
        * scale;
      
    }
      }
    } 
  }
  if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
 
    for (PetscInt k=lzs; k<lze; k++) {
      for (PetscInt i=lxs; i<lxe; i++) {
    PetscInt j=ys;
    
    PetscReal dpdc = 0.;
    PetscReal dpdz = 0.;
    if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
      if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
        dpdc = (p[k][j][i  ] + p[k][j+1][i  ] -
            p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
      }
    }
    else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
      if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
        dpdc = (p[k][j][i  ] + p[k][j+1][i  ] -
            p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
      }
    }
    else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
      if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
        dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
            p[k][j][i  ] - p[k][j+1][i  ]) * 0.5;
      }
    }
    else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
      if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
        dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
            p[k][j][i  ] - p[k][j+1][i  ]) * 0.5;
      }
    }
    else {
      dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
          p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
    }
    
    PetscReal dpde = p[k][j+1][i] - p[k][j][i];
    
    if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
      if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
        dpdz = (p[k  ][j][i] + p[k  ][j+1][i] -
            p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
      }
    }
    else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
      if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
        dpdz = (p[k  ][j][i] + p[k  ][j+1][i] -
            p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
      }
    }
    else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
      if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
        dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
            p[k  ][j][i] - p[k  ][j+1][i]) * 0.5;
      }
    }
    else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
      if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
        dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
            p[k  ][j][i] - p[k  ][j+1][i]) * 0.5;
      }
    }
    else {
      dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
          p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
    }
        
    if (!(nvert[k][j][i] + nvert[k][j+1][i])) {
      ucont[k][j][i].y -=
        (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
             jcsi[k][j][i].y * jeta[k][j][i].y +
             jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
         dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
             jeta[k][j][i].y * jeta[k][j][i].y +
             jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
         dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
             jzet[k][j][i].y * jeta[k][j][i].y +
             jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i])
        * scale;
    }
      }
    }
  }
 
  if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
    for (PetscInt j=lys; j<lye; j++) {
      for (PetscInt i=lxs; i<lxe; i++) {
    
    PetscInt k=zs;
    PetscReal dpdc = 0.;
    PetscReal dpde = 0.;
 
    if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
      if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
        dpdc = (p[k][j][i  ] + p[k+1][j][i  ] -
            p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
      }
    }
    else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
      if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
        dpdc = (p[k][j][i  ] + p[k+1][j][i  ] -
            p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
      }
    }
    else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
      if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
        dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
            p[k][j][i  ] - p[k+1][j][i  ]) * 0.5;
      }
    }
    else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
      if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
        dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
            p[k][j][i  ] - p[k+1][j][i  ]) * 0.5;
      }
    }
    else {
      dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
          p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
    }
    
    if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
      if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
        dpde = (p[k][j  ][i] + p[k+1][j  ][i] -
            p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
      }
    }
    else  if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
      if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
        dpde = (p[k][j  ][i] + p[k+1][j  ][i] -
            p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
      }
    }
    else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
      if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
        dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
            p[k][j  ][i] - p[k+1][j  ][i]) * 0.5;
      }
    }
    else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
      if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
        dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
            p[k][j  ][i] - p[k+1][j  ][i]) * 0.5;
      }
    }
    else {
      dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
          p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
    }
    
    PetscReal dpdz = p[k+1][j][i] - p[k][j][i];
    
    if (!(nvert[k][j][i] + nvert[k+1][j][i])) {
      
      ucont[k][j][i].z -=
       (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
                   kcsi[k][j][i].y * kzet[k][j][i].y +
                   kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
               dpde * (keta[k][j][i].x * kzet[k][j][i].x +
                   keta[k][j][i].y * kzet[k][j][i].y +
                   keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
               dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
                   kzet[k][j][i].y * kzet[k][j][i].y +
                   kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i])
                        * scale;
      
    }
      }
    }
  }
  
  // Corrects Flux Profile for Driven Flows if applicable.
  CorrectChannelFluxProfile(user);
  
  //================================================================================
  // Section 3: Finalization and Cleanup
  //================================================================================
 
  // --- Restore access to all PETSc vector arrays ---
  DMDAVecRestoreArray(fda, user->Ucont, &ucont);
  // DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
  //DMDAVecRestoreArray(da, user->lAj, &aj);
  DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
  DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
  DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
  DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
  DMDAVecRestoreArray(da, user->lPhi, &p);
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
 
  // --- Update ghost cells for the newly corrected velocity field ---
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating ghost cells for corrected velocity.\n");
  const FieldId staggered_fields[] = {FIELD_ID_UCONT};
  ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
  // --- Convert velocity to Cartesian and update ghost nodes ---
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Converting velocity to Cartesian and finalizing ghost nodes.\n");
  ierr = Contra2Cart(user); CHKERRQ(ierr);
  ierr = FinalizePostProjectionCellFields(user); CHKERRQ(ierr);
  //GhostNodeVelocity(user);
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting Projection step.\n");
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "UpdatePressure"
/**
 * @brief Implementation of \ref UpdatePressure().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see UpdatePressure()
 */
PetscErrorCode UpdatePressure(UserCtx *user)
{
  PetscErrorCode ierr;
 
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering UpdatePressure.\n");
 
  //================================================================================
  // Section 1: Initialization and Data Acquisition
  //================================================================================
  DM da = user->da;
  DMDALocalInfo info = user->info;
 
  // Local grid extents for the main update loop
  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;
 
  // --- Get direct pointer access to PETSc vector data for performance ---
  PetscReal ***p, ***phi;
  DMDAVecGetArray(da, user->P, &p);
  DMDAVecGetArray(da, user->Phi, &phi);
 
  //================================================================================
  // Section 2: Core Pressure Update
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Performing core pressure update (P_new = P_old + Phi).\n");
  for (PetscInt k = zs; k < ze; k++) {
    for (PetscInt j = ys; j < ye; j++) {
      for (PetscInt i = xs; i < xe; i++) {
        // This is the fundamental pressure update in a projection method.
        p[k][j][i] += phi[k][j][i];
      }
    }
  }
 
  // Restore arrays now that the core computation is done.
  DMDAVecRestoreArray(da, user->Phi, &phi);
  DMDAVecRestoreArray(da, user->P, &p);
 
  
  //================================================================================
  // Section 3: Handle Periodic Boundary Condition Synchronization
  //================================================================================
  const FieldId periodic_fields[] = {FIELD_ID_P, FIELD_ID_PHI};
  ierr = SynchronizePeriodicCellFields(user, 2, periodic_fields); CHKERRQ(ierr);
  
  //================================================================================
  // Section 4: Final Cleanup (pointers already restored)
  //================================================================================
 
  ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
  ierr = UpdateLocalGhosts(user, FIELD_ID_PHI); CHKERRQ(ierr);
  
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting UpdatePressure.\n");
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "PoissonNullSpaceFunction"
/**
 * @brief Implementation of \ref PoissonNullSpaceFunction().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see PoissonNullSpaceFunction()
 */
 
PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp,Vec X, void *ctx)
{
  PetscErrorCode ierr;
  UserCtx *user = (UserCtx*)ctx;
  (void)nullsp;
 
  DM da = user->da;
 
  DMDALocalInfo info = user->info;
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
  PetscInt  lxs, lxe, lys, lye, lzs, lze;
 
  PetscReal ***x, ***nvert;
  PetscInt  i, j, k;
 
/*   /\* First remove a constant from the Vec field X *\/ */
 
 
  /* Then apply boundary conditions */
  DMDAVecGetArray(da, X, &x);
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
  PetscReal lsum, sum;
  PetscReal  lnum, num;
 
  if (user->multinullspace) {
    LOG_ALLOW(GLOBAL, LOG_INFO, "Poisson solve is using the configured multi-nullspace.\n");
  }
  if (!user->multinullspace) {
    lsum = 0;
    lnum = 0;
    for (k=lzs; k<lze; k++) {
      for (j=lys; j<lye; j++) {
    for (i=lxs; i<lxe; i++) {
      if (nvert[k][j][i] < 0.1) {
        lsum += x[k][j][i];
        lnum ++;
      }
    }
      }
    }
 
    ierr = MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
   /*  PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
/*     PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
    sum = sum / (-1.0 * num);
 
    for (k=lzs; k<lze; k++) {
      for (j=lys; j<lye; j++) {
    for (i=lxs; i<lxe; i++) {
      if (nvert[k][j][i] < 0.1) {
        x[k][j][i] +=sum;
      }
    }
      }
    }
  }
  else {
    lsum = 0;
    lnum = 0;
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    for (k=lzs; k<lze; k++) {
      if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
        lsum += x[k][j][i];
        lnum ++;
      }
    }
      }
    }
    ierr = MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
   /*  PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
/*     PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
    sum /= -num;
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    for (k=lzs; k<lze; k++) {
      if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
        x[k][j][i] += sum;
      }
    }
      }
    }
 
    lsum = 0;
    lnum = 0;
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    for (k=lzs; k<lze; k++) {
      if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
        lsum += x[k][j][i];
        lnum ++;
      }
    }
      }
    }
    ierr = MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
   /*  PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
/*     PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
    sum /= -num;
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    for (k=lzs; k<lze; k++) {
      if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
        x[k][j][i] += sum;
      }
    }
      }
    }
               
  } //if multinullspace
  if (zs == 0) {
    k = 0;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  if (ze == mz) {
    k = mz-1;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  if (ys == 0) {
    j = 0;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  if (ye == my) {
    j = my-1;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  if (xs == 0) {
    i = 0;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  if (xe == mx) {
    i = mx-1;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    x[k][j][i] = 0.;
      }
    }
  }
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (nvert[k][j][i] > 0.1)
      x[k][j][i] = 0.;
      }
    }
  }
  DMDAVecRestoreArray(da, X, &x);
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
 
  return 0;
}
 
/**
 * @brief Implementation of \ref MyInterpolation().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see MyInterpolation()
 */
 
PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
{
  UserCtx *user;
 
  MatShellGetContext(A, (void**)&user);
 
  
  
  DM    da = user->da;
 
  DM    da_c = *user->da_c;
 
  DMDALocalInfo info = user->info;
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
  PetscInt  lxs, lxe, lys, lye, lzs, lze;
  
  PetscReal ***f, ***x, ***nvert, ***nvert_c;
  PetscInt i, j, k, ic, jc, kc, ia, ja, ka;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
 
  DMDAVecGetArray(da,   F, &f);
 
 
  Vec lX;
  DMCreateLocalVector(da_c, &lX);
 
  DMGlobalToLocalBegin(da_c, X, INSERT_VALUES, lX);
  DMGlobalToLocalEnd(da_c, X, INSERT_VALUES, lX);  
  DMDAVecGetArray(da_c, lX, &x);
 
  DMDAVecGetArray(da, user->lNvert, &nvert);
  DMDAVecGetArray(da_c, *(user->lNvert_c), &nvert_c);
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
 
    GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user);
 
      f[k][j][i] = (x[kc   ][jc   ][ic   ] * 9 +
            x[kc   ][jc+ja][ic   ] * 3 +
            x[kc   ][jc   ][ic+ia] * 3 +
            x[kc   ][jc+ja][ic+ia]) * 3./64. +
        (x[kc+ka][jc   ][ic   ] * 9 +
         x[kc+ka][jc+ja][ic   ] * 3 +
         x[kc+ka][jc   ][ic+ia] * 3 +
         x[kc+ka][jc+ja][ic+ia]) /64.;
      }
    }
  }
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
 
    if (i==0) {
      f[k][j][i] = 0.;//-f[k][j][i+1];
    }
    else if (i==mx-1) {
      f[k][j][i] = 0.;//-f[k][j][i-1];
    }
    else if (j==0) {
      f[k][j][i] = 0.;//-f[k][j+1][i];
    }
    else if (j==my-1) {
      f[k][j][i] = 0.;//-f[k][j-1][i];
    }
    else if (k==0) {
      f[k][j][i] = 0.;//-f[k+1][j][i];
    }
    else if (k==mz-1) {
      f[k][j][i] = 0.;//-f[k-1][j][i];
    }
    if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
 
      }
    }
  }
 
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
  DMDAVecRestoreArray(da_c, *(user->lNvert_c), &nvert_c);
 
  DMDAVecRestoreArray(da_c, lX, &x);
 
  VecDestroy(&lX);
  DMDAVecRestoreArray(da,   F,  &f);
 
 
 
  return 0;
 
}
 
/**
 * @brief Restrict residuals while accounting for solid-cell occupancy in the stencil.
 */
static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
{
  UserCtx *user;
  MatShellGetContext(A, (void**)&user);
  
  DM    da = user->da;
  DM    da_f = *user->da_f;
 
  DMDALocalInfo info;
  DMDAGetLocalInfo(da, &info);
  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;
  
  PetscReal ***f, ***x, ***nvert;
  PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
 
  DMDAVecGetArray(da,   F, &f);
 
  Vec lX;
  DMCreateLocalVector(da_f, &lX);
  DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
  DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);  
  DMDAVecGetArray(da_f, lX, &x);
 
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  PetscReal ***nvert_f;
  DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
 
  if ((user->isc)) ia = 0;
  else ia = 1;
 
  if ((user->jsc)) ja = 0;
  else ja = 1;
 
  if ((user->ksc)) ka = 0;
  else ka = 1;
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
        // --- CORRECTED LOGIC ---
        // First, check if the current point is a boundary point.
        // If it is, it does not contribute to the coarse grid residual.
        if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1 || nvert[k][j][i] > 0.1) {
            f[k][j][i] = 0.0;
        } 
        // Only if it's a true interior fluid point, perform the restriction.
        else {
            GridRestriction(i, j, k, &ih, &jh, &kh, user);
            f[k][j][i] = 0.125 *
              (x[kh   ][jh   ][ih   ] * PetscMax(0., 1 - nvert_f[kh   ][jh   ][ih   ]) +
               x[kh   ][jh   ][ih-ia] * PetscMax(0., 1 - nvert_f[kh   ][jh   ][ih-ia]) +
               x[kh   ][jh-ja][ih   ] * PetscMax(0., 1 - nvert_f[kh   ][jh-ja][ih   ]) +
               x[kh-ka][jh   ][ih   ] * PetscMax(0., 1 - nvert_f[kh-ka][jh   ][ih   ]) +
               x[kh   ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh   ][jh-ja][ih-ia]) +
               x[kh-ka][jh-ja][ih   ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih   ]) +
               x[kh-ka][jh   ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh   ][ih-ia]) +
               x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
        }
      }
    }
  }
 
  DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
  DMDAVecRestoreArray(da_f, lX, &x);
  VecDestroy(&lX);
  DMDAVecRestoreArray(da,   F,  &f);
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
 
  return 0;
}
 
/**
 * @brief Implementation of \ref MyRestriction().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see MyRestriction()
 */
 
PetscErrorCode MyRestriction(Mat A, Vec X, Vec F)
{
  UserCtx *user;
 
  MatShellGetContext(A, (void**)&user);
 
  
  DM    da = user->da;
 
  DM    da_f = *user->da_f;
 
  DMDALocalInfo info;
  DMDAGetLocalInfo(da, &info);
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
  //  PetscInt  lxs, lxe, lys, lye, lzs, lze;
  
  PetscReal ***f, ***x, ***nvert;
  PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
 
  DMDAVecGetArray(da,   F, &f);
 
  Vec lX;
 
  DMCreateLocalVector(da_f, &lX);
  DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
  DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);  
  DMDAVecGetArray(da_f, lX, &x);
 
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  PetscReal ***nvert_f;
  DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
 
  if ((user->isc)) ia = 0;
  else ia = 1;
 
  if ((user->jsc)) ja = 0;
  else ja = 1;
 
  if ((user->ksc)) ka = 0;
  else ka = 1;
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (k==0) {
      f[k][j][i] = 0.;
    }
    else if (k==mz-1) {
      f[k][j][i] = 0.;
    }
    else if (j==0) {
      f[k][j][i] = 0.;
    }
    else if (j==my-1) {
      f[k][j][i] = 0.;
    }
    else if (i==0) {
      f[k][j][i] = 0.;
    }
    else if (i==mx-1) {
      f[k][j][i] = 0.;
    }
    else {
      GridRestriction(i, j, k, &ih, &jh, &kh, user);
      f[k][j][i] = 0.125 *
        (x[kh   ][jh   ][ih   ] * PetscMax(0., 1 - nvert_f[kh   ][jh   ][ih   ]) +
         x[kh   ][jh   ][ih-ia] * PetscMax(0., 1 - nvert_f[kh   ][jh   ][ih-ia]) +
         x[kh   ][jh-ja][ih   ] * PetscMax(0., 1 - nvert_f[kh   ][jh-ja][ih   ]) +
         x[kh-ka][jh   ][ih   ] * PetscMax(0., 1 - nvert_f[kh-ka][jh   ][ih   ]) +
         x[kh   ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh   ][jh-ja][ih-ia]) +
         x[kh-ka][jh-ja][ih   ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih   ]) +
         x[kh-ka][jh   ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh   ][ih-ia]) +
         x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
 
 
 
      if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
    }
      }
    }
  }
 
 
  DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
 
  DMDAVecRestoreArray(da_f, lX, &x);
  VecDestroy(&lX);
 
  DMDAVecRestoreArray(da,   F,  &f);
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
 
 
  return 0;
}
 
#undef __FUNCT__
#define __FUNCT__ "PoissonLHSNew"
 
/**
 * @brief Internal helper implementation: `PoissonLHSNew()`.
 * @details Local to this translation unit.
 */
PetscErrorCode PoissonLHSNew(UserCtx *user)
{
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonLHSNew to assemble Laplacian matrix.\n");
  PetscErrorCode ierr;
  //================================================================================
  // Section 1: Initialization and Data Acquisition
  //================================================================================
 
  
  // --- Get simulation and grid context ---
  DM da = user->da, fda = user->fda;
  DMDALocalInfo info = user->info;
  PetscInt IM = user->IM, JM = user->JM, KM = user->KM;
  PetscInt i,j,k;
 
  // --- Grid dimensions ---
  PetscInt mx = info.mx, my = info.my, mz = info.mz;
  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 gxs = info.gxs, gxe = gxs + info.gxm;
  PetscInt gys = info.gys, gye = gys + info.gym;
  PetscInt gzs = info.gzs, gze = gzs + info.gzm;
 
  // --- Define constants for clarity ---
  const PetscReal IBM_FLUID_THRESHOLD = 0.1;
 
  // --- Allocate the LHS matrix A on the first call ---
  if (!user->assignedA) {
    LOG_ALLOW(GLOBAL, LOG_INFO, "First call: Creating LHS matrix 'A' with 19-point stencil preallocation.\n");
    PetscInt N = mx * my * mz; // Total size
    PetscInt M;                // Local size
    VecGetLocalSize(user->Phi, &M);
    // Create a sparse AIJ matrix, preallocating for 19 non-zeros per row (d=diagonal, o=off-diagonal)
    MatCreateAIJ(PETSC_COMM_WORLD, M, M, N, N, 19, PETSC_NULLPTR, 19, PETSC_NULLPTR, &(user->A));
    user->assignedA = PETSC_TRUE;
  }
 
  // Zero out matrix entries from the previous solve
  MatZeroEntries(user->A);
 
  // --- Get direct pointer access to grid metric data ---
  Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
  PetscReal ***aj, ***iaj, ***jaj, ***kaj, ***nvert;
  DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
  DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
  DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
  DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
  DMDAVecGetArray(da, user->lAj, &aj); DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  // --- Create temporary vectors for the metric tensor components G_ij ---
  Vec G11, G12, G13, G21, G22, G23, G31, G32, G33;
  PetscReal ***g11, ***g12, ***g13, ***g21, ***g22, ***g23, ***g31, ***g32, ***g33;
  VecDuplicate(user->lAj, &G11); VecDuplicate(user->lAj, &G12); VecDuplicate(user->lAj, &G13);
  VecDuplicate(user->lAj, &G21); VecDuplicate(user->lAj, &G22); VecDuplicate(user->lAj, &G23);
  VecDuplicate(user->lAj, &G31); VecDuplicate(user->lAj, &G32); VecDuplicate(user->lAj, &G33);
  DMDAVecGetArray(da, G11, &g11); DMDAVecGetArray(da, G12, &g12); DMDAVecGetArray(da, G13, &g13);
  DMDAVecGetArray(da, G21, &g21); DMDAVecGetArray(da, G22, &g22); DMDAVecGetArray(da, G23, &g23);
  DMDAVecGetArray(da, G31, &g31); DMDAVecGetArray(da, G32, &g32); DMDAVecGetArray(da, G33, &g33);
 
  //================================================================================
  // Section 2: Pre-compute Metric Tensor Coefficients (g_ij)
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pre-computing metric tensor components (g_ij).\n");
  for (k = gzs; k < gze; k++) {
    for (j = gys; j < gye; j++) {
      for (i = gxs; i < gxe; i++) {
        // These coefficients represent the dot products of the grid's contravariant base vectors,
        // scaled by face area. They are the core of the Laplacian operator on a curvilinear grid.
        if(i>-1 && j>-1 && k>-1 && i<IM+1 && j<JM+1 && k<KM+1){
            g11[k][j][i] = (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
            g12[k][j][i] = (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
            g13[k][j][i] = (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
            g21[k][j][i] = (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
            g22[k][j][i] = (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
            g23[k][j][i] = (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
            g31[k][j][i] = (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
            g32[k][j][i] = (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
            g33[k][j][i] = (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
        }
      }
    }
  }
 
  //================================================================================
  // Section 3: Assemble the LHS Matrix A
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling the LHS matrix A using a 19-point stencil.\n");
 
  // --- Define domain boundaries for stencil logic, accounting for periodic BCs ---
  PetscInt x_str, x_end, y_str, y_end, z_str, z_end;
  if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) { x_end = mx - 1; x_str = 0; }
  else { x_end = mx - 2; x_str = 1; }
  if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) { y_end = my - 1; y_str = 0; }
  else { y_end = my - 2; y_str = 1; }
  if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) { z_end = mz - 1; z_str = 0; }
  else { z_end = mz - 2; z_str = 1; }
 
  // --- Main assembly loop over all local grid points ---
  for (k = zs; k < ze; k++) {
    for (j = ys; j < ye; j++) {
      for (i = xs; i < xe; i++) {
        PetscScalar vol[19]; // Holds the 19 stencil coefficient values for the current row
        PetscInt idx[19];    // Holds the 19 global column indices for the current row
        PetscInt row = Gidx(i, j, k, user); // Global index for the current row
 
        // --- Handle Domain Boundary and Immersed Solid Points ---
        // For these points, we don't solve the Poisson equation. We set an identity
        // row (A_ii = 1) to effectively fix the pressure value (usually to 0).
        if (i == 0 || i == mx - 1 || j == 0 || j == my - 1 || k == 0 || k == mz - 1 || nvert[k][j][i] > IBM_FLUID_THRESHOLD) {
          vol[CP] = 1.0;
          idx[CP] = row;
          MatSetValues(user->A, 1, &row, 1, &idx[CP], &vol[CP], INSERT_VALUES);
        }
        // --- Handle Fluid Points ---
        else {
          for (PetscInt m = 0; m < 19; m++) {
            vol[m] = 0.0;
          }
 
          /************************************************************************
           * EAST FACE CONTRIBUTION (between i and i+1)
           ************************************************************************/
          if (nvert[k][j][i + 1] < IBM_FLUID_THRESHOLD && i != x_end) { // East neighbor is fluid
            // Primary derivative term: d/d_csi (g11 * dP/d_csi)
            vol[CP] -= g11[k][j][i];
            vol[EP] += g11[k][j][i];
 
            // Cross-derivative term: d/d_csi (g12 * dP/d_eta).
            // This requires an average of dP/d_eta. If a neighbor is solid, the stencil
            // dynamically switches to a one-sided difference to avoid using solid points.
          if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
          vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
          vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
        }
          }
          else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
          vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
          vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
        }
          }
          else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
          vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
          vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
        }
          }
          else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
          vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
          vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
        }
          }
          else { // Centered difference
        vol[NP] += g12[k][j][i] * 0.25; vol[NE] += g12[k][j][i] * 0.25;
        vol[SP] -= g12[k][j][i] * 0.25; vol[SE] -= g12[k][j][i] * 0.25;
          }
            
            // Cross-derivative term: d/d_csi (g13 * dP/d_zet)
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
          vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
          vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
        }
          }
          else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
          vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
          vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
        }
          }
          else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
          vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
          vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
        }
          }
          else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
          vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
          vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
        }
          }
          else { // Centered difference
        vol[TP] += g13[k][j][i] * 0.25; vol[TE] += g13[k][j][i] * 0.25;
        vol[BP] -= g13[k][j][i] * 0.25; vol[BE] -= g13[k][j][i] * 0.25;
          }
          }
 
          /************************************************************************
           * WEST FACE CONTRIBUTION (between i-1 and i)
           ************************************************************************/
          if (nvert[k][j][i-1] < IBM_FLUID_THRESHOLD && i != x_str) {
          vol[CP] -= g11[k][j][i-1];
          vol[WP] += g11[k][j][i-1];
 
          if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
          vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
          vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
        }
          }
          else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1) {
          vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
          vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
        }
          }
          else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
          vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
          vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
        }
          }
          else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
          vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
          vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
        }
          }
          else {
        vol[NP] -= g12[k][j][i-1] * 0.25; vol[NW] -= g12[k][j][i-1] * 0.25;
        vol[SP] += g12[k][j][i-1] * 0.25; vol[SW] += g12[k][j][i-1] * 0.25;
          }
 
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
          vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
          vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
        }
          }
          else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1) {
          vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
          vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
        }
          }
          else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
          vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
          vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
        }
          }
          else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
          vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
          vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
        }
          }
          else {
        vol[TP] -= g13[k][j][i-1] * 0.25; vol[TW] -= g13[k][j][i-1] * 0.25;
        vol[BP] += g13[k][j][i-1] * 0.25; vol[BW] += g13[k][j][i-1] * 0.25;
          }
          }
 
          /************************************************************************
           * NORTH FACE CONTRIBUTION (between j and j+1)
           ************************************************************************/
          if (nvert[k][j+1][i] < IBM_FLUID_THRESHOLD && j != y_end) {
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
          vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
          vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
        }
          }
          else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
          vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
          vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
        }
          }
          else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
          vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
          vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
        }
          }
          else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC &&  nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
          vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
          vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
        }
          }
          else {
        vol[EP] += g21[k][j][i] * 0.25; vol[NE] += g21[k][j][i] * 0.25;
        vol[WP] -= g21[k][j][i] * 0.25; vol[NW] -= g21[k][j][i] * 0.25;
          }
 
          vol[CP] -= g22[k][j][i];
          vol[NP] += g22[k][j][i];
 
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
          vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
          vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
        }
          }
          else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
          vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
          vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
        }
          }
          else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
          vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
          vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
        }
          }
          else if ((k == 1 || k==mz-2 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
          vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
          vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
        }
          }
          else {
        vol[TP] += g23[k][j][i] * 0.25; vol[TN] += g23[k][j][i] * 0.25;
        vol[BP] -= g23[k][j][i] * 0.25; vol[BN] -= g23[k][j][i] * 0.25;
          }
          }
 
          /************************************************************************
           * SOUTH FACE CONTRIBUTION (between j-1 and j)
           ************************************************************************/
          if (nvert[k][j-1][i] < IBM_FLUID_THRESHOLD && j != y_str) {
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
          vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
          vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
        }
          }
          else  if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1) {
          vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
          vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
        }
          }
          else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
          vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
          vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
        }
          }
          else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
          vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
          vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
        }
          }
          else {
        vol[EP] -= g21[k][j-1][i] * 0.25; vol[SE] -= g21[k][j-1][i] * 0.25;
        vol[WP] += g21[k][j-1][i] * 0.25; vol[SW] += g21[k][j-1][i] * 0.25;
          }
 
          vol[CP] -= g22[k][j-1][i];
          vol[SP] += g22[k][j-1][i];
 
          if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
          vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
          vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
        }
          }
          else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 ) {
          vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
          vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
        }
          }
          else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
          vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
          vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
        }
          }
          else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
        if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
          vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
          vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
        }
          }
          else {
        vol[TP] -= g23[k][j-1][i] * 0.25; vol[TS] -= g23[k][j-1][i] * 0.25;
        vol[BP] += g23[k][j-1][i] * 0.25; vol[BS] += g23[k][j-1][i] * 0.25;
          }
          }
 
          /************************************************************************
           * TOP FACE CONTRIBUTION (between k and k+1)
           ************************************************************************/
          if (nvert[k+1][j][i] < IBM_FLUID_THRESHOLD && k != z_end) {
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
          vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
          vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
        }
          }
          else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
          vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
          vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
        }
          }
          else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
          vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
          vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
        }
          }
          else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
          vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
          vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
        }
          }
          else {
        vol[EP] += g31[k][j][i] * 0.25; vol[TE] += g31[k][j][i] * 0.25;
        vol[WP] -= g31[k][j][i] * 0.25; vol[TW] -= g31[k][j][i] * 0.25;
          }
 
          if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
          vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
          vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
        }
          }
          else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
          vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
          vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
        }
          }
          else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
          vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
          vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
        }
          }
          else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
          vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
          vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
        }
          }
          else {
        vol[NP] += g32[k][j][i] * 0.25; vol[TN] += g32[k][j][i] * 0.25;
        vol[SP] -= g32[k][j][i] * 0.25; vol[TS] -= g32[k][j][i] * 0.25;
          }
 
          vol[CP] -= g33[k][j][i];
          vol[TP] += g33[k][j][i];
          }
 
          /************************************************************************
           * BOTTOM FACE CONTRIBUTION (between k-1 and k)
           ************************************************************************/
          if (nvert[k-1][j][i] < IBM_FLUID_THRESHOLD && k != z_str) {
          if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
          vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
          vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
        }
          }
          else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
        if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1) {
          vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
          vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
        }
          }
          else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
          vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
          vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
        }
          }
          else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
        if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
          vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
          vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
        }
          }
          else {
        vol[EP] -= g31[k-1][j][i] * 0.25; vol[BE] -= g31[k-1][j][i] * 0.25;
        vol[WP] += g31[k-1][j][i] * 0.25; vol[BW] += g31[k-1][j][i] * 0.25;
          }
 
          if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
          vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
          vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
        }
          }
          else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
        if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1) {
          vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
          vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
        }
          }
          else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
          vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
          vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
        }
          }
          else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
        if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
          vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
          vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
        }
          }
          else {
        vol[NP] -= g32[k-1][j][i] * 0.25; vol[BN] -= g32[k-1][j][i] * 0.25;
        vol[SP] += g32[k-1][j][i] * 0.25; vol[BS] += g32[k-1][j][i] * 0.25;
          }
 
          vol[CP] -= g33[k-1][j][i];
          vol[BP] += g33[k-1][j][i];
          }
 
          // --- Final scaling and insertion into the matrix ---
 
          // Scale all stencil coefficients by the negative cell volume (-aj).
          for (PetscInt m = 0; m < 19; m++) {
            vol[m] *= -aj[k][j][i];
          }
 
          // Set the global column indices for the 19 stencil points, handling periodic BCs.
          idx[CP] = Gidx(i, j, k, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[EP] = Gidx(1, j, k, user); else idx[EP] = Gidx(i+1, j, k, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[WP] = Gidx(mx-2, j, k, user); else idx[WP] = Gidx(i-1, j, k, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NP] = Gidx(i, 1, k, user); else idx[NP] = Gidx(i, j+1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SP] = Gidx(i, my-2, k, user); else idx[SP] = Gidx(i, j-1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TP] = Gidx(i, j, 1, user); else idx[TP] = Gidx(i, j, k+1, user);
          if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BP] = Gidx(i, j, mz-2, user); else idx[BP] = Gidx(i, j, k-1, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==mx-2 && j==my-2) idx[NE] = Gidx(1, 1, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[NE] = Gidx(1, j+1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NE] = Gidx(i+1, 1, k, user); else idx[NE] = Gidx(i+1, j+1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==mx-2 && j==1) idx[SE] = Gidx(1, my-2, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[SE] = Gidx(1, j-1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SE] = Gidx(i+1, my-2, k, user); else idx[SE] = Gidx(i+1, j-1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==1 && j==my-2) idx[NW] = Gidx(mx-2, 1, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[NW] = Gidx(mx-2, j+1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[NW] = Gidx(i-1, 1, k, user); else idx[NW] = Gidx(i-1, j+1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && i==1 && j==1) idx[SW] = Gidx(mx-2, my-2, k, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[SW] = Gidx(mx-2, j-1, k, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[SW] = Gidx(i-1, my-2, k, user); else idx[SW] = Gidx(i-1, j-1, k, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==my-2 && k==mz-2) idx[TN] = Gidx(i, 1, 1, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[TN] = Gidx(i, 1, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TN] = Gidx(i, j+1, 1, user); else idx[TN] = Gidx(i, j+1, k+1, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==my-2 && k==1) idx[BN] = Gidx(i, 1, mz-2, user); else if(user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2) idx[BN] = Gidx(i, 1, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BN] = Gidx(i, j+1, mz-2, user); else idx[BN] = Gidx(i, j+1, k-1, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==1 && k==mz-2) idx[TS] = Gidx(i, my-2, 1, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[TS] = Gidx(i, my-2, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TS] = Gidx(i, j-1, 1, user); else idx[TS] = Gidx(i, j-1, k+1, user);
          if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && j==1 && k==1) idx[BS] = Gidx(i, my-2, mz-2, user); else if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==1) idx[BS] = Gidx(i, my-2, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BS] = Gidx(i, j-1, mz-2, user); else idx[BS] = Gidx(i, j-1, k-1, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==mx-2 && k==mz-2) idx[TE] = Gidx(1, j, 1, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[TE] = Gidx(1, j, k+1, user); else if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TE] = Gidx(i+1, j, 1, user); else idx[TE] = Gidx(i+1, j, k+1, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==mx-2 && k==1) idx[BE] = Gidx(1, j, mz-2, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2) idx[BE] = Gidx(1, j, k-1, user); else if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BE] = Gidx(i+1, j, mz-2, user); else idx[BE] = Gidx(i+1, j, k-1, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==1 && k==mz-2) idx[TW] = Gidx(mx-2, j, 1, user); else if(user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[TW] = Gidx(mx-2, j, k+1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2) idx[TW] = Gidx(i-1, j, 1, user); else idx[TW] = Gidx(i-1, j, k+1, user);
          if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && i==1 && k==1) idx[BW] = Gidx(mx-2, j, mz-2, user); else if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==1) idx[BW] = Gidx(mx-2, j, k-1, user); else if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==1) idx[BW] = Gidx(i-1, j, mz-2, user); else idx[BW] = Gidx(i-1, j, k-1, user);
 
          // Insert the computed row into the matrix A.
          MatSetValues(user->A, 1, &row, 19, idx, vol, INSERT_VALUES);
        }
      }
    }
  }
 
  //================================================================================
  // Section 4: Finalize Matrix and Cleanup
  //================================================================================
  
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing matrix assembly.\n");
  MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY);
 
  PetscReal max_A;
 
  ierr = MatNorm(user->A,NORM_INFINITY,&max_A);CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG," Max value in A matrix for level %d =  %le.\n",user->thislevel,max_A);
 
  // if (get_log_level() >= LOG_DEBUG) {
  //  ierr = MatView(user->A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  // }
  
  // --- Restore access to all PETSc vectors and destroy temporary ones ---
  DMDAVecRestoreArray(da, G11, &g11); DMDAVecRestoreArray(da, G12, &g12); DMDAVecRestoreArray(da, G13, &g13);
  DMDAVecRestoreArray(da, G21, &g21); DMDAVecRestoreArray(da, G22, &g22); DMDAVecRestoreArray(da, G23, &g23);
  DMDAVecRestoreArray(da, G31, &g31); DMDAVecRestoreArray(da, G32, &g32); DMDAVecRestoreArray(da, G33, &g33);
  
  VecDestroy(&G11); VecDestroy(&G12); VecDestroy(&G13);
  VecDestroy(&G21); VecDestroy(&G22); VecDestroy(&G23);
  VecDestroy(&G31); VecDestroy(&G32); VecDestroy(&G33);
 
  DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
  DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
  DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
  DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
  DMDAVecRestoreArray(da, user->lAj, &aj); DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
  
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting PoissonLHSNew.\n");
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "PoissonRHS"
 
/**
 * @brief Implementation of \ref PoissonRHS().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see PoissonRHS()
 */
 
PetscErrorCode PoissonRHS(UserCtx *user, Vec B)
{
  PetscErrorCode ierr;
  DMDALocalInfo info = user->info;
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt      i, j, k;
  PetscReal ***nvert, ***aj, ***rb, dt = user->simCtx->dt;
  struct Components{
    PetscReal x;
    PetscReal y;
    PetscReal z;
  } *** ucont;
 
  PROFILE_FUNCTION_BEGIN;
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonRHS to compute pressure equation RHS.\n");
 
  DMDAVecGetArray(user->da, B, &rb);
  DMDAVecGetArray(user->fda, user->lUcont, &ucont);
  DMDAVecGetArray(user->da, user->lNvert, &nvert);
  DMDAVecGetArray(user->da, user->lAj, &aj);
 
 
   LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing RHS values for each cell.\n");
  
  for (k=zs; k<ze; k++) { 
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
 
    if (i==0 || i==mx-1 || j==0 || j==my-1 ||  k==0 || k==mz-1) {
      rb[k][j][i] = 0.;
    }
    else if (nvert[k][j][i] > 0.1) {
      rb[k][j][i] = 0;
    }
    else {
      rb[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
              ucont[k][j][i].y - ucont[k][j-1][i].y +
              ucont[k][j][i].z - ucont[k-1][j][i].z) / dt
        * aj[k][j][i] / 1.0 * COEF_TIME_ACCURACY;     // user->simCtx->st replaced by 1.0.
     
    }
      }
    }
  }
 
 
  // --- Check the solvability condition for the Poisson equation ---
  // The global sum of the RHS (proportional to the total divergence) must be zero.
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying solvability condition (sum of RHS terms).\n");
  PetscReal lsum=0., sum=0.;
 
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    
    lsum += rb[k][j][i] / aj[k][j][i]* dt/COEF_TIME_ACCURACY;
 
      }
    }
  }
  
  ierr = MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
  LOG_ALLOW(GLOBAL, LOG_INFO, "Global Sum of RHS (Divergence Check): %le\n", sum);
 
  user->simCtx->summationRHS = sum;
    
  DMDAVecRestoreArray(user->fda, user->lUcont, &ucont);
  DMDAVecRestoreArray(user->da, user->lNvert, &nvert);
  DMDAVecRestoreArray(user->da, user->lAj, &aj);
  DMDAVecRestoreArray(user->da, B, &rb);
 
 
  PROFILE_FUNCTION_END;
  return 0;
}
 
/**
 * @brief Implementation of \ref VolumeFlux_rev().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see VolumeFlux_rev()
 */
 
PetscErrorCode VolumeFlux_rev(UserCtx *user, PetscReal *ibm_Flux, 
                  PetscReal *ibm_Area, PetscInt flg)
{
  PetscErrorCode ierr;
 
  DM    da = user->da, fda = user->fda;
 
  DMDALocalInfo info = user->info;
 
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt i, j, k;
  PetscInt  lxs, lys, lzs, lxe, lye, lze;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
  PetscReal ***nvert, ibmval=1.5;
  Cmpnts ***ucor, ***csi, ***eta, ***zet;
  DMDAVecGetArray(fda, user->Ucont, &ucor);
  DMDAVecGetArray(fda, user->lCsi, &csi);
  DMDAVecGetArray(fda, user->lEta, &eta);
  DMDAVecGetArray(fda, user->lZet, &zet);
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  PetscReal libm_Flux, libm_area;
  libm_Flux = 0;
  libm_area = 0;
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
        libm_Flux += ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
        libm_Flux += ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
        libm_Flux += ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
    if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
      if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
        libm_Flux -= ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] < 0.1 && j < my-2) {
        libm_Flux -= ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
        libm_Flux -= ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
      }
    }
  }
 
  ierr = MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
  ierr = MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
 /*  PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
/*   PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "IBM flux correction: flux=%le, area=%le\n", *ibm_Flux, *ibm_Area);
 
  PetscReal correction;
 
  if (*ibm_Area > 1.e-15) {
    if (flg)
      correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
    else
      correction = *ibm_Flux / *ibm_Area;
  }
  else {
    correction = 0;
  }
 
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
        ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
                    correction;
              
      }
      if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
        ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x + 
                    eta[k][j][i].y * eta[k][j][i].y +
                    eta[k][j][i].z * eta[k][j][i].z) *
                    correction;
      }
      if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
        ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
                    zet[k][j][i].y * zet[k][j][i].y +
                    zet[k][j][i].z * zet[k][j][i].z) *
                    correction;
      }
    }
 
    if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
      if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
        ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
                    correction;
              
      }
      if (nvert[k][j+1][i] < 0.1 && j < my-2) {
        ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                    eta[k][j][i].y * eta[k][j][i].y +
                    eta[k][j][i].z * eta[k][j][i].z) *
                    correction;
      }
      if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
        ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                    zet[k][j][i].y * zet[k][j][i].y +
                    zet[k][j][i].z * zet[k][j][i].z) *
                    correction;
      }
    }
 
      }
    }
  }
  
 
 
  libm_Flux = 0;
  libm_area = 0;
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
        libm_Flux += ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
        libm_Flux += ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
        libm_Flux += ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
    if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
      if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
        libm_Flux -= ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] < 0.1 && j < my-2) {
        libm_Flux -= ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
        libm_Flux -= ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
      }
    }
  }
 
  ierr = MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
  ierr = MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
 /*  PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
/*   PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "IBM flux measurement: flux=%le, area=%le\n", *ibm_Flux, *ibm_Area);
 
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
  DMDAVecRestoreArray(fda, user->lCsi, &csi);
  DMDAVecRestoreArray(fda, user->lEta, &eta);
  DMDAVecRestoreArray(fda, user->lZet, &zet);
  DMDAVecRestoreArray(fda, user->Ucont, &ucor);
 
  const FieldId staggered_fields[] = {FIELD_ID_UCONT};
  ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
  return 0;
}
 
 
/**
 * @brief Implementation of \ref VolumeFlux().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see VolumeFlux()
 */
 
PetscErrorCode VolumeFlux(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
{
  PetscErrorCode ierr;
  // --- CONTEXT ACQUISITION BLOCK ---
  // Get the master simulation context from the UserCtx.
  SimCtx *simCtx = user->simCtx;
 
  // Create local variables to mirror the legacy globals for minimal code changes.
  const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
  // --- END CONTEXT ACQUISITION BLOCK ---
 
  DM    da = user->da, fda = user->fda;
 
  DMDALocalInfo info = user->info;
 
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt i, j, k,ibi;
  PetscInt  lxs, lys, lzs, lxe, lye, lze;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
  
  PetscReal epsilon=1.e-8;
  PetscReal ***nvert, ibmval=1.9999;
  
  struct Components {
    PetscReal x;
    PetscReal y;
    PetscReal z;
  }***ucor, ***csi, ***eta, ***zet;
 
 
  PetscInt xend=mx-2 ,yend=my-2,zend=mz-2;
 
  if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) xend=mx-1;
  if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) yend=my-1;
  if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) zend=mz-1;
 
  DMDAVecGetArray(fda, user->Ucont, &ucor);
  DMDAVecGetArray(fda, user->lCsi, &csi);
  DMDAVecGetArray(fda, user->lEta, &eta);
  DMDAVecGetArray(fda, user->lZet, &zet);
  DMDAVecGetArray(da, user->lNvert, &nvert);
 
  PetscReal libm_Flux, libm_area, libm_Flux_abs=0., ibm_Flux_abs;
  libm_Flux = 0;
  libm_area = 0;
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering VolumeFlux to enforce no-penetration condition.\n");
 
  //Mohsen March 2017
  PetscReal *lIB_Flux = NULL, *lIB_area = NULL, *IB_Flux = NULL, *IB_Area = NULL;
  if (NumberOfBodies > 1) { 
  
    lIB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
    lIB_area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
    IB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
    IB_Area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
  
 
    for (ibi=0; ibi<NumberOfBodies; ibi++) {
      lIB_Flux[ibi]=0.0;
      lIB_area[ibi]=0.0;
      IB_Flux[ibi]=0.0;
      IB_Area[ibi]=0.0;
    }
  }
 
 
  //================================================================================
  // PASS 1: Calculate Uncorrected Flux and Area
  // This pass measures the total fluid "leakage" across the immersed boundary.
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 1: Measuring uncorrected flux and area.\n");
  
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
        
        if (fabs(ucor[k][j][i].x)>epsilon) {
          libm_Flux += ucor[k][j][i].x;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                csi[k][j][i].y * csi[k][j][i].y +
                                csi[k][j][i].z * csi[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].x);
          
          libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                csi[k][j][i].y * csi[k][j][i].y +
                csi[k][j][i].z * csi[k][j][i].z);
          
          if (NumberOfBodies > 1) {
        
        ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
        lIB_Flux[ibi] += ucor[k][j][i].x;
        lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                      csi[k][j][i].y * csi[k][j][i].y +
                      csi[k][j][i].z * csi[k][j][i].z);
          }
        } else
          ucor[k][j][i].x=0.;
        
      }
      if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
        
        if (fabs(ucor[k][j][i].y)>epsilon) {
          libm_Flux += ucor[k][j][i].y;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].y)/sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                eta[k][j][i].y * eta[k][j][i].y +
                                eta[k][j][i].z * eta[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].y);
          libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                eta[k][j][i].y * eta[k][j][i].y +
                eta[k][j][i].z * eta[k][j][i].z);
          if (NumberOfBodies > 1) {
        
        ibi=(int)((nvert[k][j+1][i]-1.0)*1001); 
 
        lIB_Flux[ibi] += ucor[k][j][i].y;
        lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                      eta[k][j][i].y * eta[k][j][i].y +
                      eta[k][j][i].z * eta[k][j][i].z);
          }
        } else
          ucor[k][j][i].y=0.;
      }
      if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
        
        if (fabs(ucor[k][j][i].z)>epsilon) {
          libm_Flux += ucor[k][j][i].z;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                zet[k][j][i].y * zet[k][j][i].y +
                                zet[k][j][i].z * zet[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].z);
          libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                zet[k][j][i].y * zet[k][j][i].y +
                zet[k][j][i].z * zet[k][j][i].z);
          
          if (NumberOfBodies > 1) {
        
        ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
        lIB_Flux[ibi] += ucor[k][j][i].z;
        lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                      zet[k][j][i].y * zet[k][j][i].y +
                      zet[k][j][i].z * zet[k][j][i].z);
          }
        }else
          ucor[k][j][i].z=0.;
      }
    }
 
    if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
      
      if (nvert[k][j][i+1] < 0.1 && i < xend) {
        if (fabs(ucor[k][j][i].x)>epsilon) {
          libm_Flux -= ucor[k][j][i].x;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                csi[k][j][i].y * csi[k][j][i].y +
                                csi[k][j][i].z * csi[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].x);
          libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                csi[k][j][i].y * csi[k][j][i].y +
                csi[k][j][i].z * csi[k][j][i].z);
          if (NumberOfBodies > 1) {
        
        ibi=(int)((nvert[k][j][i]-1.0)*1001);
        lIB_Flux[ibi] -= ucor[k][j][i].x;
        lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                      csi[k][j][i].y * csi[k][j][i].y +
                      csi[k][j][i].z * csi[k][j][i].z);
          }
              
        }else
          ucor[k][j][i].x=0.;
      }
      if (nvert[k][j+1][i] < 0.1 && j < yend) {
        if (fabs(ucor[k][j][i].y)>epsilon) {
          libm_Flux -= ucor[k][j][i].y;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].y)/ sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                 eta[k][j][i].y * eta[k][j][i].y +
                                 eta[k][j][i].z * eta[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].y);
          libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                eta[k][j][i].y * eta[k][j][i].y +
                eta[k][j][i].z * eta[k][j][i].z);
          if (NumberOfBodies > 1) {
 
        ibi=(int)((nvert[k][j][i]-1.0)*1001);
        lIB_Flux[ibi] -= ucor[k][j][i].y;
        lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                      eta[k][j][i].y * eta[k][j][i].y +
                      eta[k][j][i].z * eta[k][j][i].z);
          }
        }else
          ucor[k][j][i].y=0.;
      }
      if (nvert[k+1][j][i] < 0.1 && k < zend) {
        if (fabs(ucor[k][j][i].z)>epsilon) {
          libm_Flux -= ucor[k][j][i].z;
          if (flg==3)
        libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                zet[k][j][i].y * zet[k][j][i].y +
                                zet[k][j][i].z * zet[k][j][i].z);
          else
        libm_Flux_abs += fabs(ucor[k][j][i].z);
          libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                zet[k][j][i].y * zet[k][j][i].y +
                zet[k][j][i].z * zet[k][j][i].z);
          if (NumberOfBodies > 1) {
        
        ibi=(int)((nvert[k][j][i]-1.0)*1001);
        lIB_Flux[ibi] -= ucor[k][j][i].z;
        lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                      zet[k][j][i].y * zet[k][j][i].y +
                      zet[k][j][i].z * zet[k][j][i].z);
          }
        }else
          ucor[k][j][i].z=0.;
      }
    }
    
      }
    }
  }
  
  ierr = MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
  ierr = MPI_Allreduce(&libm_Flux_abs, &ibm_Flux_abs,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
  ierr = MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
  if (NumberOfBodies > 1) { 
    ierr = MPI_Allreduce(lIB_Flux,IB_Flux,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
    ierr = MPI_Allreduce(lIB_area,IB_Area,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRMPI(ierr);
  }
 
  PetscReal correction;
 
  PetscReal *Correction = NULL;
  if (NumberOfBodies > 1) {
      Correction=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
      for (ibi=0; ibi<NumberOfBodies; ibi++) Correction[ibi]=0.0;
  }
 
  if (*ibm_Area > 1.e-15) {
    if (flg>1) 
      correction = (*ibm_Flux + user->FluxIntpSum)/ ibm_Flux_abs;
    else if (flg)
      correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
    else
      correction = *ibm_Flux / *ibm_Area;
    if (NumberOfBodies > 1) 
      for (ibi=0; ibi<NumberOfBodies; ibi++) if (IB_Area[ibi]>1.e-15) Correction[ibi] = IB_Flux[ibi] / IB_Area[ibi];
  }
  else {
    correction = 0;
  }
  // --- Log the uncorrected results and calculated correction ---
  LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Uncorrected Flux: %g, Area: %g, Correction: %g\n", *ibm_Flux, *ibm_Area, correction);
  if  (NumberOfBodies>1){
    for (ibi=0; ibi<NumberOfBodies; ibi++)   LOG_ALLOW(GLOBAL, LOG_INFO, "  [Body %d] Uncorrected Flux: %g, Area: %g, Correction: %g\n", ibi, IB_Flux[ibi], IB_Area[ibi], Correction[ibi]);
  }
 
  //================================================================================
  // PASS 2: Apply Correction to Velocity Field
  // This pass modifies the velocity at the boundary to enforce zero net flux.
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 2: Applying velocity corrections at the boundary.\n");
  
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] <ibmval && i < xend) {
        if (fabs(ucor[k][j][i].x)>epsilon){
          if (flg==3) 
        ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x)/
          sqrt(csi[k][j][i].x * csi[k][j][i].x +
               csi[k][j][i].y * csi[k][j][i].y +
               csi[k][j][i].z * csi[k][j][i].z);
          else if (flg==2) 
        ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x);
          else if (NumberOfBodies > 1) {
        ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
        ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
          Correction[ibi];
          }
          else
        ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
          correction;
        }
      }
      if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
        if (fabs(ucor[k][j][i].y)>epsilon) {
          if (flg==3) 
        ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y)/
          sqrt(eta[k][j][i].x * eta[k][j][i].x + 
               eta[k][j][i].y * eta[k][j][i].y +
               eta[k][j][i].z * eta[k][j][i].z);
          else if (flg==2) 
        ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y);
          else if (NumberOfBodies > 1) {
        ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
        ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
                    eta[k][j][i].y * eta[k][j][i].y +
                    eta[k][j][i].z * eta[k][j][i].z) *
          Correction[ibi];
          }
          else
        ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x + 
                    eta[k][j][i].y * eta[k][j][i].y +
                    eta[k][j][i].z * eta[k][j][i].z) *
          correction;
        }
      }
      if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
        if (fabs(ucor[k][j][i].z)>epsilon) {
          if (flg==3) 
        ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z)/
          sqrt(zet[k][j][i].x * zet[k][j][i].x +
               zet[k][j][i].y * zet[k][j][i].y +
               zet[k][j][i].z * zet[k][j][i].z);
          else if (flg==2) 
        ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z);
          else if (NumberOfBodies > 1) {
        ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
        ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
                    zet[k][j][i].y * zet[k][j][i].y +
                    zet[k][j][i].z * zet[k][j][i].z) *
          Correction[ibi];
          }
          else
        ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
                    zet[k][j][i].y * zet[k][j][i].y +
                    zet[k][j][i].z * zet[k][j][i].z) *
          correction;
        }
      }
    }
 
    if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
      if (nvert[k][j][i+1] < 0.1 && i < xend) {
        if (fabs(ucor[k][j][i].x)>epsilon) {
          if (flg==3) 
        ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x)/
          sqrt(csi[k][j][i].x * csi[k][j][i].x +
               csi[k][j][i].y * csi[k][j][i].y +
               csi[k][j][i].z * csi[k][j][i].z);
          else if (flg==2) 
        ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x);
          else if (NumberOfBodies > 1) {
        ibi=(int)((nvert[k][j][i]-1.0)*1001);
        ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
          Correction[ibi];
          }
          else
        ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                    csi[k][j][i].y * csi[k][j][i].y +
                    csi[k][j][i].z * csi[k][j][i].z) *
          correction;
        }
      }
      if (nvert[k][j+1][i] < 0.1 && j < yend) {
        if (fabs(ucor[k][j][i].y)>epsilon) {
        if (flg==3) 
          ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y)/
        sqrt(eta[k][j][i].x * eta[k][j][i].x +
             eta[k][j][i].y * eta[k][j][i].y +
             eta[k][j][i].z * eta[k][j][i].z);
        else if (flg==2) 
          ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y);
        else if (NumberOfBodies > 1) {
          ibi=(int)((nvert[k][j][i]-1.0)*1001);
          ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                      eta[k][j][i].y * eta[k][j][i].y +
                      eta[k][j][i].z * eta[k][j][i].z) *
        Correction[ibi];
        }
        else
        ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                    eta[k][j][i].y * eta[k][j][i].y +
                    eta[k][j][i].z * eta[k][j][i].z) *
                    correction;
        }
      }
      if (nvert[k+1][j][i] < 0.1 && k < zend) {
        if (fabs(ucor[k][j][i].z)>epsilon) {
        if (flg==3) 
          ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z)/
        sqrt(zet[k][j][i].x * zet[k][j][i].x +
             zet[k][j][i].y * zet[k][j][i].y +
             zet[k][j][i].z * zet[k][j][i].z);
        else if (flg==2) 
          ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z);
        else if (NumberOfBodies > 1) {
          ibi=(int)((nvert[k][j][i]-1.0)*1001);
          ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                      zet[k][j][i].y * zet[k][j][i].y +
                      zet[k][j][i].z * zet[k][j][i].z) *
        Correction[ibi];
        }
        else
          ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                      zet[k][j][i].y * zet[k][j][i].y +
                      zet[k][j][i].z * zet[k][j][i].z) *
        correction;
        }
      }
    }
    
      }
    }
  }
  
  //================================================================================
  // PASS 3: Verification
  // This optional pass recalculates the flux to confirm the correction was successful.
  //================================================================================
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 3: Verifying corrected flux.\n");
  
  libm_Flux = 0;
  libm_area = 0;
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
        libm_Flux += ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
        libm_Flux += ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
        libm_Flux += ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
    if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
      if (nvert[k][j][i+1] < 0.1 && i < xend) {
        libm_Flux -= ucor[k][j][i].x;
        libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
              csi[k][j][i].y * csi[k][j][i].y +
              csi[k][j][i].z * csi[k][j][i].z);
              
      }
      if (nvert[k][j+1][i] < 0.1 && j < yend) {
        libm_Flux -= ucor[k][j][i].y;
        libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
              eta[k][j][i].y * eta[k][j][i].y +
              eta[k][j][i].z * eta[k][j][i].z);
      }
      if (nvert[k+1][j][i] < 0.1 && k < zend) {
        libm_Flux -= ucor[k][j][i].z;
        libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
              zet[k][j][i].y * zet[k][j][i].y +
              zet[k][j][i].z * zet[k][j][i].z);
      }
    }
 
      }
    }
  }
 
  ierr = MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
  ierr = MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
 /*  PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
/*   PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
    LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Corrected (Verified) Flux: %g, Area: %g\n", *ibm_Flux, *ibm_Area);
 
 
  if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
    if (xe==mx){
      i=mx-2;
      for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      // if(j>0 && k>0 && j<user->JM && k<user->KM){
        if ((nvert[k][j][i]>ibmval && nvert[k][j][i+1]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j][i+1]>ibmval)) ucor[k][j][i].x=0.0;
        
        // }
    }
      }
    }
  }
 
  if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
    if (ye==my){
      j=my-2;
      for (k=lzs; k<lze; k++) {
    for (i=lxs; i<lxe; i++) {
      // if(i>0 && k>0 && i<user->IM && k<user->KM){
        if ((nvert[k][j][i]>ibmval && nvert[k][j+1][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j+1][i]>ibmval)) ucor[k][j][i].y=0.0;
        // }
    }
      }
    }
  }
 
  if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
    if (ze==mz){
      k=mz-2;
      for (j=lys; j<lye; j++) {
    for (i=lxs; i<lxe; i++) {
      // if(i>0 && j>0 && i<user->IM && j<user->JM){
        if ((nvert[k][j][i]>ibmval && nvert[k+1][j][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k+1][j][i]>ibmval)) ucor[k][j][i].z=0.0;
        // }
    }
      }
    }
  }
 
 
  DMDAVecRestoreArray(da, user->lNvert, &nvert);
  DMDAVecRestoreArray(fda, user->lCsi, &csi);
  DMDAVecRestoreArray(fda, user->lEta, &eta);
  DMDAVecRestoreArray(fda, user->lZet, &zet);
  DMDAVecRestoreArray(fda, user->Ucont, &ucor);
 
  const FieldId staggered_fields[] = {FIELD_ID_UCONT};
  ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
  if (NumberOfBodies > 1) {
    free(lIB_Flux);
    free(lIB_area);
    free(IB_Flux);
    free(IB_Area);
    free(Correction);
  }
 
 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting VolumeFlux.\n");
  
  return 0;
}
 
/**
 * @brief Report whether a coarse-grid cell is completely blocked by solid fine-grid cells.
 */
static PetscErrorCode FullyBlocked(UserCtx *user)
{
  PetscErrorCode ierr;
  DM da = user->da;
  Vec nNvert;
  DMDALocalInfo info = user->info;
 
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt i, j, k;
 
  PetscInt *KSKE = user->KSKE;
  PetscReal ***nvert;
  PetscBool *Blocked;
 
  DMDACreateNaturalVector(da, &nNvert);
  DMDAGlobalToNaturalBegin(da, user->Nvert, INSERT_VALUES, nNvert);
  DMDAGlobalToNaturalEnd(da, user->Nvert, INSERT_VALUES, nNvert);
 
  VecScatter ctx;
  Vec Zvert;
  VecScatterCreateToZero(nNvert, &ctx, &Zvert);
 
  VecScatterBegin(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
  VecScatterEnd(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
 
  VecScatterDestroy(&ctx);
  VecDestroy(&nNvert);
 
  PetscInt rank;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
 
  if (!rank) {
 
    VecGetArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
    PetscMalloc(mx*my*sizeof(PetscBool), &Blocked);
    for (j=1; j<my-1; j++) {
      for (i=1; i<mx-1; i++) {
    Blocked[j*mx+i] = PETSC_FALSE;
    for (k=0; k<mz; k++) {
      if (nvert[k][j][i] > 0.1) {
        if (!Blocked[j*mx+i]) {
          KSKE[2*(j*mx+i)] = k;
          Blocked[j*mx+i] = PETSC_TRUE;
        }
        else {
          KSKE[2*(j*mx+i)] = PetscMin(KSKE[2*(j*mx+i)], k);
        }
      }
    }
      }
    }
 
 
    user->multinullspace = PETSC_TRUE;
    for (j=1; j<my-1; j++) {
      for (i=1; i<mx-1; i++) {
    if (!Blocked[j*mx+i]) {
      user->multinullspace = PETSC_FALSE;
      break;
    }
      }
    }
    PetscFree(Blocked);
    VecRestoreArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
    ierr = MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    if (user->multinullspace) {
      ierr = MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    }
  }
  else {
    ierr = MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    if (user->multinullspace) {
      ierr = MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRMPI(ierr);
    }
  }
 
 
 
  VecDestroy(&Zvert);
  return 0;
}
 
/**
 * @brief Restrict solid-volume fractions from fine cells to one coarse cell.
 */
static PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
{
  PetscErrorCode ierr;
  //  DA        da = user_c->da, fda = user_c->fda;
 
 
 
  DMDALocalInfo info = user_c->info;
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt i,j,k;
  PetscInt ih, jh, kh, ia, ja, ka;
  PetscInt  lxs, lxe, lys, lye, lzs, lze;
 
  PetscReal ***nvert, ***nvert_h;
 
  DMDAVecGetArray(user_h->da, user_h->lNvert, &nvert_h);
  DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert);
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
  if ((user_c->isc)) ia = 0;
  else ia = 1;
 
  if ((user_c->jsc)) ja = 0;
  else ja = 1;
 
  if ((user_c->ksc)) ka = 0;
  else ka = 1;
 
  VecSet(user_c->Nvert, 0.);
  if (user_c->thislevel > 0) {
    for (k=lzs; k<lze; k++) {
      for (j=lys; j<lye; j++) {
    for (i=lxs; i<lxe; i++) {
      GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
      if (nvert_h[kh   ][jh   ][ih   ] *
          nvert_h[kh   ][jh   ][ih-ia] *
          nvert_h[kh   ][jh-ja][ih   ] *
          nvert_h[kh-ka][jh   ][ih   ] *
          nvert_h[kh   ][jh-ja][ih-ia] *
          nvert_h[kh-ka][jh   ][ih-ia] *
          nvert_h[kh-ka][jh-ja][ih   ] *
          nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
        nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
      }
    }
      }
    }
  }
  else {
    for (k=lzs; k<lze; k++) {
      for (j=lys; j<lye; j++) {
    for (i=lxs; i<lxe; i++) {
      GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
      if (nvert_h[kh   ][jh   ][ih   ] *
          nvert_h[kh   ][jh   ][ih-ia] *
          nvert_h[kh   ][jh-ja][ih   ] *
          nvert_h[kh-ka][jh   ][ih   ] *
          nvert_h[kh   ][jh-ja][ih-ia] *
          nvert_h[kh-ka][jh   ][ih-ia] *
          nvert_h[kh-ka][jh-ja][ih   ] *
          nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
        nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
      }
    }
      }
    }
  }
  DMDAVecRestoreArray(user_h->da, user_h->lNvert, &nvert_h);
  DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert);
 
  ierr = UpdateLocalGhosts(user_c, FIELD_ID_NVERT); CHKERRQ(ierr);
  //Mohsen Dec 2015
  DMDAVecGetArray(user_c->da, user_c->lNvert, &nvert);
  DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert_h);
 
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    if (nvert_h[k][j][i] < 0.1) {
      if (nvert[k][j][i+1] + nvert[k][j][i-1] > 1.1 &&
          nvert[k][j+1][i] + nvert[k][j-1][i] > 1.1 &&
          nvert[k+1][j][i] + nvert[k-1][j][i] > 1.1) {
        nvert_h[k][j][i] = 1.;
      }
    }
      }
    }
  }
 
  DMDAVecRestoreArray(user_c->da, user_c->lNvert, &nvert);
  DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert_h);
  ierr = UpdateLocalGhosts(user_c, FIELD_ID_NVERT); CHKERRQ(ierr);
 /*  DMLocalToGlobalBegin(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
/*   DMLocalToGlobalEnd(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
  return 0;
}
 
#undef __FUNCT__
#define __FUNCT__ "PoissonSolver_MG"
/**
 * @brief Implementation of \ref PoissonSolver_MG().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/poisson.h`.
 * @see PoissonSolver_MG()
 */
 
PetscErrorCode PoissonSolver_MG(UserMG *usermg)
{
    // --- CONTEXT ACQUISITION BLOCK ---
    // Get the master simulation context from the first block's UserCtx on the finest level.
    // This provides access to all former global variables.
    SimCtx *simCtx = usermg->mgctx[0].user[0].simCtx;
 
    // Create local variables to mirror the legacy globals for minimal code changes.
    const PetscInt block_number = simCtx->block_number;
    const PetscInt immersed = simCtx->immersed;
    const PetscInt MHV = simCtx->MHV;
    const PetscInt LV = simCtx->LV;
    PetscMPIInt rank = simCtx->rank;
    // --- END CONTEXT ACQUISITION BLOCK ---
 
    PetscErrorCode ierr;
    PetscInt l, bi;
    MGCtx *mgctx = usermg->mgctx;
    KSP mgksp, subksp;
    PC mgpc, subpc;
    UserCtx *user;
 
    PetscFunctionBeginUser; // Moved to after variable declarations
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting Multigrid Poisson Solve...\n");
 
    for (bi = 0; bi < block_number; bi++) {
        
        // ====================================================================
        //   SECTION: Immersed Boundary Specific Setup (Conditional)
        // ====================================================================
        if (immersed) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Performing IBM pre-solve setup (Nvert restriction, etc.).\n", bi);
            for (l = usermg->mglevels - 1; l > 0; l--) {
                mgctx[l].user[bi].multinullspace = PETSC_FALSE;
                MyNvertRestriction(&mgctx[l].user[bi], &mgctx[l-1].user[bi]);
            }
            // Coarsest level check for disconnected domains due to IBM
            l = 0;
            user = mgctx[l].user;
            /* KSKE is allocated once by CreateAndInitializeAllVectors; FullyBlocked
             * rewrites every entry it reads, so no per-solve reallocation is needed. */
            FullyBlocked(&user[bi]);
        }
        
 
        l = usermg->mglevels - 1;
        user = mgctx[l].user;
        
        // We are solving the  linear system AX=B where A = Laplacian Operator Matrix; X = Unknown Phi (Pressure Correction) and B = RHS (Flux Divergence based)
 
        // --- 1. Compute RHS of the Poisson Equation ---
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Computing Poisson RHS...\n", bi);
        ierr = VecDuplicate(user[bi].P, &user[bi].B); CHKERRQ(ierr);
        
        PetscReal ibm_Flux, ibm_Area;
        PetscInt flg = immersed - 1;
 
        // Calculate volume flux source terms (often from IBM)
        VolumeFlux(&user[bi], &ibm_Flux, &ibm_Area, flg);
        if (MHV || LV) {
            flg = ((MHV > 1 || LV) && bi == 0) ? 1 : 0;
            VolumeFlux_rev(&user[bi], &ibm_Flux, &ibm_Area, flg);
        }
        // Calculate the main flux divergence term B.
        PoissonRHS(&user[bi], user[bi].B);
        
        // --- 2. Assemble LHS Matrix (Laplacian) on all MG levels ---
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Assembling Poisson LHS on all levels...\n", bi);
        for (l = usermg->mglevels - 1; l >= 0; l--) {
            user = mgctx[l].user;
        LOG_ALLOW(GLOBAL,LOG_DEBUG," Calculating LHS for Level %d.\n",l);
            PoissonLHSNew(&user[bi]);
        }
 
        // --- 3. Setup PETSc KSP and PCMG (Multigrid Preconditioner) ---
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Configuring KSP and PCMG...\n", bi);
 
    ierr = KSPCreate(PETSC_COMM_WORLD, &mgksp); CHKERRQ(ierr);
        ierr = KSPAppendOptionsPrefix(mgksp, "ps_"); CHKERRQ(ierr);
 
    // =======================================================================
        DualMonitorCtx *monctx;
        char           filen[PETSC_MAX_PATH_LEN + 128];
 
        // 1. Allocate the context and set it up.
        ierr = PetscNew(&monctx); CHKERRQ(ierr);
 
    monctx->step = simCtx->step;
    monctx->block_id = bi;
    monctx->file_handle = NULL;
 
    // Only rank 0 handles the file.
        if (!rank) {
          ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Poisson_Solver_Convergence_History_Block_%d.log", simCtx->log_dir, bi); CHKERRQ(ierr);
      // On the very first step of a fresh run, TRUNCATE the file.
      // In continue mode, always APPEND to preserve existing data.
      if (simCtx->step == simCtx->StartStep + 1 && !simCtx->continueMode) {
        monctx->file_handle = fopen(filen, "w");
      } else { // For all subsequent steps (or continue mode), APPEND.
        monctx->file_handle = fopen(filen, "a");
      }
 
            if (monctx->file_handle) {
                if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
                    PetscFPrintf(PETSC_COMM_SELF, monctx->file_handle,
                                 "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
                }
                PetscFPrintf(PETSC_COMM_SELF, monctx->file_handle, "--- Convergence for Timestep %d, Block %d ---\n", (int)simCtx->step, bi);
            } else {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open KSP monitor log file: %s", filen);
            }
        }
 
        monctx->log_to_console = simCtx->ps_ksp_pic_monitor_true_residual;
 
        ierr = KSPMonitorSet(mgksp, DualKSPMonitor, monctx, DualMonitorDestroy); CHKERRQ(ierr);
        // =======================================================================
    
        ierr = KSPGetPC(mgksp, &mgpc); CHKERRQ(ierr);
        ierr = PCSetType(mgpc, PCMG); CHKERRQ(ierr);
 
        ierr = PCMGSetLevels(mgpc, usermg->mglevels, PETSC_NULLPTR); CHKERRQ(ierr);
        ierr = PCMGSetCycleType(mgpc, PC_MG_CYCLE_V); CHKERRQ(ierr);
        ierr = PCMGSetType(mgpc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr);
        if (simCtx->mg_preItr != simCtx->mg_poItr) {
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "PETSc PCMG exposes one smoother count in this build; using max(pre_sweeps=%d, post_sweeps=%d).\n",
                      simCtx->mg_preItr, simCtx->mg_poItr);
        }
        PetscInt mg_smooths = simCtx->mg_preItr > simCtx->mg_poItr ? simCtx->mg_preItr : simCtx->mg_poItr;
        ierr = PCMGSetNumberSmooth(mgpc, mg_smooths); CHKERRQ(ierr);
 
        // --- 4. Define Restriction and Interpolation Operators for MG ---
        for (l = usermg->mglevels - 1; l > 0; l--) {
 
         // Get stable pointers directly from the main mgctx array.
         // These pointers point to memory that will persist.
      UserCtx *fine_user_ctx   = &mgctx[l].user[bi];
      UserCtx *coarse_user_ctx = &mgctx[l-1].user[bi];
 
      // --- Configure the context pointers ---
      // The coarse UserCtx needs to know about the fine grid for restriction.
      coarse_user_ctx->da_f     = &(fine_user_ctx->da);
      coarse_user_ctx->user_f   = fine_user_ctx;
 
      // The fine UserCtx needs to know about the coarse grid for interpolation.
      fine_user_ctx->da_c       = &(coarse_user_ctx->da);
      fine_user_ctx->user_c     = coarse_user_ctx;
      fine_user_ctx->lNvert_c   = &(coarse_user_ctx->lNvert);
 
      // --- Get matrix dimensions ---
      PetscInt m_c = (coarse_user_ctx->info.xm * coarse_user_ctx->info.ym * coarse_user_ctx->info.zm);
      PetscInt m_f = (fine_user_ctx->info.xm * fine_user_ctx->info.ym * fine_user_ctx->info.zm);
      PetscInt M_c = (coarse_user_ctx->info.mx * coarse_user_ctx->info.my * coarse_user_ctx->info.mz);
      PetscInt M_f = (fine_user_ctx->info.mx * fine_user_ctx->info.my * fine_user_ctx->info.mz);
 
      LOG_ALLOW(GLOBAL,LOG_DEBUG,"level = %d; m_c = %d; m_f = %d; M_c = %d; M_f = %d.\n",l,m_c,m_f,M_c,M_f);
      // --- Create the MatShell objects ---
      // Pass the STABLE pointer coarse_user_ctx as the context for restriction.
      ierr = MatCreateShell(PETSC_COMM_WORLD, m_c, m_f, M_c, M_f, (void*)coarse_user_ctx, &fine_user_ctx->MR); CHKERRQ(ierr);
    
      // Pass the STABLE pointer fine_user_ctx as the context for interpolation.
      ierr = MatCreateShell(PETSC_COMM_WORLD, m_f, m_c, M_f, M_c, (void*)fine_user_ctx, &fine_user_ctx->MP); CHKERRQ(ierr);
    
      // --- Set the operations for the MatShells ---
      ierr = MatShellSetOperation(fine_user_ctx->MR, MATOP_MULT, (void(*)(void))RestrictResidual_SolidAware); CHKERRQ(ierr);
      ierr = MatShellSetOperation(fine_user_ctx->MP, MATOP_MULT, (void(*)(void))MyInterpolation); CHKERRQ(ierr);
    
      // --- Register the operators with PCMG ---
      ierr = PCMGSetRestriction(mgpc, l, fine_user_ctx->MR); CHKERRQ(ierr);
      ierr = PCMGSetInterpolation(mgpc, l, fine_user_ctx->MP); CHKERRQ(ierr);
      
        }
        
        // --- 5. Configure Solvers on Each MG Level ---
        for (l = usermg->mglevels - 1; l >= 0; l--) {
            user = mgctx[l].user;
            if (l > 0) { // Smoother for fine levels
                ierr = PCMGGetSmoother(mgpc, l, &subksp); CHKERRQ(ierr);
            } else { // Direct or iterative solver for the coarsest level
                ierr = PCMGGetCoarseSolve(mgpc, &subksp); CHKERRQ(ierr);
                ierr = KSPSetTolerances(subksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, 40); CHKERRQ(ierr);
            } 
            
            ierr = KSPSetOperators(subksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
            ierr = KSPGetPC(subksp, &subpc); CHKERRQ(ierr);
        ierr = PCSetType(subpc, PCBJACOBI); CHKERRQ(ierr);
        ierr = KSPSetFromOptions(subksp); CHKERRQ(ierr);
        
        PCType subpc_type;
        PetscBool is_bjacobi = PETSC_FALSE;
        ierr = PCGetType(subpc, &subpc_type); CHKERRQ(ierr);
        if (subpc_type) {
          ierr = PetscStrcmp(subpc_type, PCBJACOBI, &is_bjacobi); CHKERRQ(ierr);
        }
 
        if (is_bjacobi) {
          KSP *subsubksp;
          PC subsubpc;
          PetscInt nlocal;
 
          ierr = KSPSetUp(subksp); CHKERRQ(ierr); // Set up KSP to allow access to sub-KSPs
          ierr = PCBJacobiGetSubKSP(subpc, &nlocal, NULL, &subsubksp); CHKERRQ(ierr);
 
          for (PetscInt abi = 0; abi < nlocal; abi++) {
            ierr = KSPGetPC(subsubksp[abi], &subsubpc); CHKERRQ(ierr);
            // Add the critical shift amount for the nested block-Jacobi factor PC.
            ierr = PCFactorSetShiftAmount(subsubpc, 1.e-10); CHKERRQ(ierr);
          }
        }
        
            ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULLPTR, &user[bi].nullsp); CHKERRQ(ierr);
            ierr = MatNullSpaceSetFunction(user[bi].nullsp, PoissonNullSpaceFunction, &user[bi]); CHKERRQ(ierr);
            ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
            
            ierr = PCMGSetResidual(mgpc, l, PCMGResidualDefault, user[bi].A); CHKERRQ(ierr);
            ierr = KSPSetUp(subksp); CHKERRQ(ierr);
 
            if (l < usermg->mglevels - 1) {
                ierr = MatCreateVecs(user[bi].A, &user[bi].R, PETSC_NULLPTR); CHKERRQ(ierr);
                ierr = PCMGSetRhs(mgpc, l, user[bi].R); CHKERRQ(ierr);
            }
        }
 
        // --- 6. Set Final KSP Operators and Solve ---
        l = usermg->mglevels - 1;
        user = mgctx[l].user;
        
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Setting KSP operators and solving...\n", bi);
        ierr = KSPSetOperators(mgksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
        ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
        ierr = KSPSetFromOptions(mgksp); CHKERRQ(ierr);
        ierr = KSPSetUp(mgksp); CHKERRQ(ierr);
        ierr = KSPSolve(mgksp, user[bi].B, user[bi].Phi); CHKERRQ(ierr);
 
        // --- 7. Cleanup for this block ---
        for (l = usermg->mglevels - 1; l >= 0; l--) {
            user = mgctx[l].user;
            MatNullSpaceDestroy(&user[bi].nullsp);
            MatDestroy(&user[bi].A);
            user[bi].assignedA = PETSC_FALSE;
            if (l > 0) {
                MatDestroy(&user[bi].MR);
                MatDestroy(&user[bi].MP);
            }
            if (l < usermg->mglevels - 1) {
                VecDestroy(&user[bi].R);
            }
        }
        
        KSPDestroy(&mgksp);
        VecDestroy(&mgctx[usermg->mglevels-1].user[bi].B);
        
    } // End of loop over blocks
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Multigrid Poisson Solve complete.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}