poisson.h File Reference

#include "variables.h"
#include "Metric.h"
#include "Boundaries.h"

Include dependency graph for poisson.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
PetscErrorCode	PoissonSolver_MG (UserMG *usermg)
	Solves the pressure-Poisson equation using a geometric multigrid method.
PetscErrorCode	PoissonLHSNew (UserCtx *user)
	Assembles the Left-Hand-Side (LHS) matrix (Laplacian operator) for the Poisson equation on a single grid level.
PetscErrorCode	PoissonRHS (UserCtx *user, Vec B)
	Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermediate velocity field.
PetscErrorCode	UpdatePressure (UserCtx *user)
	Updates the pressure field P with the pressure correction Phi computed by the Poisson solver.
PetscErrorCode	CorrectChannelFluxProfile (UserCtx *user)
	Enforces a constant volumetric flux profile along the entire length of a driven periodic channel.
PetscErrorCode	Projection (UserCtx *user)
	Corrects the contravariant velocity field Ucont to be divergence-free using the gradient of the pressure correction field Phi.
PetscErrorCode	PoissonNullSpaceFunction (MatNullSpace nullsp, Vec X, void *ctx)
	The callback function for PETSc's MatNullSpace object.
PetscErrorCode	MyRestriction (Mat A, Vec X, Vec F)
	The callback function for the multigrid restriction operator (MatShell).
PetscErrorCode	MyInterpolation (Mat A, Vec X, Vec F)
	The callback function for the multigrid interpolation operator (MatShell).
PetscErrorCode	VolumeFlux (UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
	Calculates the net flux across the immersed boundary surface.
PetscErrorCode	VolumeFlux_rev (UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
	A specialized version of VolumeFlux, likely for reversed normals.

Function Documentation

PoissonSolver_MG()

PetscErrorCode PoissonSolver_MG ( UserMG *  usermg )

extern

Solves the pressure-Poisson equation using a geometric multigrid method.

This function orchestrates the entire multigrid V-cycle for the pressure correction equation. It assembles the Laplacian matrix on all grid levels, sets up the KSP solvers, smoothers, restriction/interpolation operators, and executes the solve.

Parameters

usermg

The UserMG context containing the entire multigrid hierarchy.

Returns

PetscErrorCode 0 on success.

Note

Testing status: This routine is exercised in runtime smoke, but still needs deeper direct bespoke coverage for debugging and branch isolation.

Solves the pressure-Poisson equation using a geometric multigrid method.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

PoissonSolver_MG()

Definition at line 3149 of file poisson.c.

{
    // --- 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;
            ierr = PetscMalloc1(user[bi].info.mx * user[bi].info.my * 2, &user[bi].KSKE); CHKERRQ(ierr);
            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);
            } else if (l==0 && immersed) {
                PetscFree(user[bi].KSKE);
            }
            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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PoissonLHSNew()

PetscErrorCode PoissonLHSNew ( UserCtx *  user )

extern

Assembles the Left-Hand-Side (LHS) matrix (Laplacian operator) for the Poisson equation on a single grid level.

Parameters

user	The UserCtx for the grid level on which to assemble the matrix.

Returns

PetscErrorCode 0 on success.

Note

Testing status: Direct unit coverage exists for core operator assembly, but periodic and immersed-boundary stencil branches remain thinner than the Cartesian baseline.

Assembles the Left-Hand-Side (LHS) matrix (Laplacian operator) for the Poisson equation on a single grid level.

Local to this translation unit.

Definition at line 1424 of file poisson.c.

{
  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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PoissonRHS()

PetscErrorCode PoissonRHS ( UserCtx *  user, Vec  B  )

extern

Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermediate velocity field.

Parameters

user	The UserCtx for the grid level.
B	The PETSc Vec where the RHS result will be stored.

Returns

PetscErrorCode 0 on success.

Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermediate velocity field.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

PoissonRHS()

Definition at line 2033 of file poisson.c.

{
  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;
}

Here is the caller graph for this function:

UpdatePressure()

PetscErrorCode UpdatePressure ( UserCtx *  user )

extern

Updates the pressure field P with the pressure correction Phi computed by the Poisson solver.

(P = P + Phi)

Parameters

user	The UserCtx containing the P and Phi vectors.

Returns

PetscErrorCode 0 on success.

Updates the pressure field P with the pressure correction Phi computed by the Poisson solver.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

UpdatePressure()

Definition at line 855 of file poisson.c.

{
  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 char *periodic_fields[] = {"P", "Phi"};
  ierr = SynchronizePeriodicCellFields(user, 2, periodic_fields); CHKERRQ(ierr);
  
  //================================================================================
  // Section 4: Final Cleanup (pointers already restored)
  //================================================================================
 
  ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
  ierr = UpdateLocalGhosts(user, "Phi"); CHKERRQ(ierr);
  
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting UpdatePressure.\n");
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CorrectChannelFluxProfile()

PetscErrorCode CorrectChannelFluxProfile

(

UserCtx *

user

)

Enforces a constant volumetric flux profile along the entire length of a driven periodic channel.

This function is a "hard" corrector, called at the end of the projection step. The projection ensures the velocity field is divergence-free (3D continuity), but this function enforces a stricter 1D continuity condition (Flux(plane) = constant) required for physically realistic, fully-developed periodic channel/pipe flow.

The process is as follows:

1. Introspects the boundary condition handlers to detect if a DRIVEN_ flow is active and in which direction ('X', 'Y', or 'Z'). If none is found, it exits.
2. Measures the current volumetric flux through every single cross-sectional plane in the driven direction.
3. For each plane, it calculates the velocity correction required to make its flux match the global targetVolumetricFlux (which was set by the controller).
4. It applies this spatially-uniform (but plane-dependent) velocity correction directly to the ucont field, ensuring Flux(plane) = TargetFlux for all planes.

Parameters

user	The UserCtx containing the simulation state for a single block.

Returns

PetscErrorCode 0 on success.

Enforces a constant volumetric flux profile along the entire length of a driven periodic channel.

Local to this translation unit.

Definition at line 107 of file poisson.c.

{
    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);
}

Here is the caller graph for this function:

Projection()

PetscErrorCode Projection ( UserCtx *  user )

extern

Corrects the contravariant velocity field Ucont to be divergence-free using the gradient of the pressure correction field Phi.

Parameters

user	The UserCtx containing the Ucont and Phi vectors.

Returns

PetscErrorCode 0 on success.

Note

Testing status: Direct unit coverage exists for basic projection invariants; periodic and immersed-boundary correction branches remain part of the next-gap backlog.

Corrects the contravariant velocity field Ucont to be divergence-free using the gradient of the pressure correction field Phi.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

Projection()

Definition at line 328 of file poisson.c.

{
  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 char *staggered_fields[] = {"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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PoissonNullSpaceFunction()

PetscErrorCode PoissonNullSpaceFunction ( MatNullSpace  nullsp, Vec  X, void *  ctx  )

extern

The callback function for PETSc's MatNullSpace object.

This function removes the null space from the Poisson solution vector by ensuring the average pressure is zero, which is necessary for problems with pure Neumann boundary conditions.

Parameters

nullsp	The MatNullSpace context.
X	The vector to be corrected.
ctx	A void pointer to the UserCtx.

Returns

PetscErrorCode 0 on success.

The callback function for PETSc's MatNullSpace object.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

PoissonNullSpaceFunction()

Definition at line 923 of file poisson.c.

{
  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) PetscPrintf(PETSC_COMM_WORLD, "MultiNullSpace!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\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;
}

Here is the caller graph for this function:

MyRestriction()

PetscErrorCode MyRestriction ( Mat  A, Vec  X, Vec  F  )

extern

The callback function for the multigrid restriction operator (MatShell).

Defines the fine-to-coarse grid transfer for the Poisson residual.

Parameters

A	The shell matrix context.
X	The fine-grid source vector.
F	The coarse-grid destination vector.

Returns

PetscErrorCode 0 on success.

The callback function for the multigrid restriction operator (MatShell).

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

MyRestriction()

Definition at line 1318 of file poisson.c.

{
  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;
}

Here is the call graph for this function:

MyInterpolation()

PetscErrorCode MyInterpolation ( Mat  A, Vec  X, Vec  F  )

extern

The callback function for the multigrid interpolation operator (MatShell).

Defines the coarse-to-fine grid transfer for the pressure correction.

Parameters

A	The shell matrix context.
X	The coarse-grid source vector.
F	The fine-grid destination vector.

Returns

PetscErrorCode 0 on success.

The callback function for the multigrid interpolation operator (MatShell).

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

MyInterpolation()

Definition at line 1125 of file poisson.c.

{
  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;
 
}

Here is the caller graph for this function:

VolumeFlux()

PetscErrorCode VolumeFlux ( UserCtx *  user, PetscReal *  ibm_Flux, PetscReal *  ibm_Area, PetscInt  flg  )

extern

Calculates the net flux across the immersed boundary surface.

Parameters

user	The UserCtx for the grid level.
ibm_Flux	(Output) The calculated net flux.
ibm_Area	(Output) The total surface area of the IB.
flg	A flag controlling the correction behavior.

Returns

PetscErrorCode 0 on success.

Calculates the net flux across the immersed boundary surface.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

VolumeFlux()

Definition at line 2364 of file poisson.c.

{
  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 char *staggered_fields[] = {"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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

VolumeFlux_rev()

PetscErrorCode VolumeFlux_rev ( UserCtx *  user, PetscReal *  ibm_Flux, PetscReal *  ibm_Area, PetscInt  flg  )

extern

A specialized version of VolumeFlux, likely for reversed normals.

Parameters

user	The UserCtx for the grid level.
ibm_Flux	(Output) The calculated net flux.
ibm_Area	(Output) The total surface area of the IB.
flg	A flag controlling the correction behavior.

Returns

PetscErrorCode 0 on success.

A specialized version of VolumeFlux, likely for reversed normals.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/poisson.h.

See also

VolumeFlux_rev()

Definition at line 2122 of file poisson.c.

{
  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); */
  PetscPrintf(PETSC_COMM_WORLD, "IBMFlux REV %le %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); */
  PetscPrintf(PETSC_COMM_WORLD, "IBMFlux22 REV %le %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 char *staggered_fields[] = {"Ucont"};
  ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function: