grid.c File Reference

#include "grid.h"
#include "logging.h"

Include dependency graph for grid.c:

Go to the source code of this file.

Macros
#define	BBOX_TOLERANCE   1e-6
#define	__FUNCT__   "ParseAndSetGridInputs"
#define	__FUNCT__   "DefineAllGridDimensions"
#define	__FUNCT__   "InitializeSingleGridDM"
#define	__FUNCT__   "InitializeAllGridDMs"
#define	__FUNCT__   "AssignAllGridCoordinates"
#define	__FUNCT__   "SetFinestLevelCoordinates"
#define	__FUNCT__   "GenerateAndSetCoordinates"
#define	__FUNCT__   "ReadAndSetCoordinates"
#define	__FUNCT__   "RestrictCoordinates"
#define	__FUNCT__   "ComputeLocalBoundingBox"
#define	__FUNCT__   "GatherAllBoundingBoxes"
#define	__FUNCT__   "BroadcastAllBoundingBoxes"
#define	__FUNCT__   "CalculateInletProperties"
#define	__FUNCT__   "CalculateOutletProperties"
#define	__FUNCT__   "CalculateFaceCenterAndArea"

Functions
static PetscErrorCode	ParseAndSetGridInputs (UserCtx *user)
	Determines the grid source and calls the appropriate parsing routine.
PetscErrorCode	DefineAllGridDimensions (SimCtx *simCtx)
	Orchestrates the parsing and setting of grid dimensions for all blocks.
static PetscErrorCode	InitializeSingleGridDM (UserCtx *user, UserCtx *coarse_user)
	Creates the DMDA objects (da and fda) for a single UserCtx.
PetscErrorCode	InitializeAllGridDMs (SimCtx *simCtx)
	Orchestrates the creation of DMDA objects for every block and multigrid level.
static PetscErrorCode	SetFinestLevelCoordinates (UserCtx *user)
	A router that populates the coordinates for a single finest-level DMDA.
static PetscErrorCode	GenerateAndSetCoordinates (UserCtx *user)
	Programmatically generates and sets grid coordinates based on user parameters.
static PetscErrorCode	ReadAndSetCoordinates (UserCtx *user, FILE *fd)
	Reads physical coordinates from a file and populates the DMDA for a specific block.
static PetscErrorCode	RestrictCoordinates (UserCtx *coarse_user, UserCtx *fine_user)
	Populates coarse grid coordinates by restricting from a fine grid.
PetscErrorCode	AssignAllGridCoordinates (SimCtx *simCtx)
	Orchestrates the assignment of physical coordinates to all DMDA objects.
static PetscReal	ComputeStretchedCoord (PetscInt i, PetscInt N, PetscReal L, PetscReal r)
	Computes a stretched coordinate along one dimension.
PetscErrorCode	ComputeLocalBoundingBox (UserCtx *user, BoundingBox *localBBox)
	Computes the local bounding box of the grid on the current process.
PetscErrorCode	GatherAllBoundingBoxes (UserCtx *user, BoundingBox **allBBoxes)
	Gathers local bounding boxes from all MPI processes to rank 0.
PetscErrorCode	BroadcastAllBoundingBoxes (UserCtx *user, BoundingBox **bboxlist)
	Broadcasts the bounding box list from rank 0 to all ranks.
PetscErrorCode	CalculateInletProperties (UserCtx *user)
	Calculates the center and area of the primary INLET face.
PetscErrorCode	CalculateOutletProperties (UserCtx *user)
	Calculates the center and area of the primary OUTLET face.
PetscErrorCode	CalculateFaceCenterAndArea (UserCtx *user, BCFace face_id, Cmpnts *face_center, PetscReal *face_area)
	Calculates the geometric center and total area of a specified boundary face.

Macro Definition Documentation

BBOX_TOLERANCE

#define BBOX_TOLERANCE   1e-6

Definition at line 6 of file grid.c.

__FUNCT__ [1/15]

#define __FUNCT__   "ParseAndSetGridInputs"

Definition at line 9 of file grid.c.

__FUNCT__ [2/15]

#define __FUNCT__   "DefineAllGridDimensions"

Definition at line 9 of file grid.c.

__FUNCT__ [3/15]

#define __FUNCT__   "InitializeSingleGridDM"

Definition at line 9 of file grid.c.

__FUNCT__ [4/15]

#define __FUNCT__   "InitializeAllGridDMs"

Definition at line 9 of file grid.c.

__FUNCT__ [5/15]

#define __FUNCT__   "AssignAllGridCoordinates"

Definition at line 9 of file grid.c.

__FUNCT__ [6/15]

#define __FUNCT__   "SetFinestLevelCoordinates"

Definition at line 9 of file grid.c.

__FUNCT__ [7/15]

#define __FUNCT__   "GenerateAndSetCoordinates"

Definition at line 9 of file grid.c.

__FUNCT__ [8/15]

#define __FUNCT__   "ReadAndSetCoordinates"

Definition at line 9 of file grid.c.

__FUNCT__ [9/15]

#define __FUNCT__   "RestrictCoordinates"

Definition at line 9 of file grid.c.

__FUNCT__ [10/15]

#define __FUNCT__   "ComputeLocalBoundingBox"

Definition at line 9 of file grid.c.

__FUNCT__ [11/15]

#define __FUNCT__   "GatherAllBoundingBoxes"

Definition at line 9 of file grid.c.

__FUNCT__ [12/15]

#define __FUNCT__   "BroadcastAllBoundingBoxes"

Definition at line 9 of file grid.c.

__FUNCT__ [13/15]

#define __FUNCT__   "CalculateInletProperties"

Definition at line 9 of file grid.c.

__FUNCT__ [14/15]

#define __FUNCT__   "CalculateOutletProperties"

Definition at line 9 of file grid.c.

__FUNCT__ [15/15]

#define __FUNCT__   "CalculateFaceCenterAndArea"

Definition at line 9 of file grid.c.

Function Documentation

ParseAndSetGridInputs()

static PetscErrorCode ParseAndSetGridInputs ( UserCtx *  user )

static

Determines the grid source and calls the appropriate parsing routine.

This function acts as a router. It checks the global simCtx->generate_grid flag (accessed via the user->simCtx back-pointer) to decide whether to call the parser for a programmatically generated grid or for a grid defined in a file.

Parameters

user	Pointer to the UserCtx for a specific block. The function will populate the geometric fields within this struct.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 22 of file grid.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx; // Get the global context via the back-pointer
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
    if(strcmp(simCtx->eulerianSource,"analytical")==0){
        ierr = SetAnalyticalGridInfo(user); CHKERRQ(ierr);
    }else{ // eulerianSource is "solve" or "load"
        if (simCtx->generate_grid) {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: Block %d is programmatically generated. Calling generation parser.\n", simCtx->rank, user->_this);
            ierr = ReadGridGenerationInputs(user); CHKERRQ(ierr);
        } else {
            LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: Block %d is file-based. Calling file parser.\n", simCtx->rank, user->_this);
            ierr = ReadGridFile(user); CHKERRQ(ierr);
        }
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DefineAllGridDimensions()

PetscErrorCode DefineAllGridDimensions

(

SimCtx *

simCtx

)

Orchestrates the parsing and setting of grid dimensions for all blocks.

This function serves as the high-level entry point for defining the geometric properties of each grid block in the simulation. It iterates through every block defined by simCtx->block_number.

For each block, it performs two key actions:

1. It explicitly sets the block's index (_this) in the corresponding UserCtx struct for the finest multigrid level. This makes the context "self-aware".
2. It calls a helper function (ParseAndSetGridInputs) to handle the detailed work of parsing options or files to populate the rest of the geometric properties for that specific block (e.g., IM, Min_X, rx).

Parameters

simCtx

The master SimCtx, which contains the number of blocks and the UserCtx hierarchy to be configured.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 68 of file grid.c.

{
    PetscErrorCode ierr;
    PetscInt       nblk = simCtx->block_number;
    UserCtx        *finest_users;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (simCtx->usermg.mglevels == 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "MG levels not set. Cannot get finest_users.");
    }
    // Get the UserCtx array for the finest grid level
    finest_users = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Defining grid dimensions for %d blocks...\n", nblk);
 
    // Loop over each block to configure its grid dimensions and geometry.
    for (PetscInt bi = 0; bi < nblk; bi++) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: --- Configuring Geometry for Block %d ---\n", simCtx->rank, bi);
 
        // Before calling any helpers, set the block index in the context.
        // This makes the UserCtx self-aware of which block it represents.
        LOG_ALLOW(GLOBAL,LOG_DEBUG,"finest_users->_this = %d, bi = %d\n",finest_users[bi]._this,bi);
        //finest_user[bi]._this = bi;
 
        // Call the helper function for this specific block. It can now derive
        // all necessary information from the UserCtx pointer it receives.
        ierr = ParseAndSetGridInputs(&finest_users[bi]); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InitializeSingleGridDM()

static PetscErrorCode InitializeSingleGridDM ( UserCtx *  user, UserCtx *  coarse_user  )

static

Creates the DMDA objects (da and fda) for a single UserCtx.

This function is a direct adaptation of the core logic in MGDACreate. It creates the scalar (da) and vector (fda) DMs for a single grid level.

If a coarse_user context is provided, it performs the critical processor alignment calculation from the legacy code. This ensures the new (fine) DM aligns with the coarse DM for multigrid efficiency. If coarse_user is NULL, it creates the DM with a default PETSc decomposition, intended for the coarsest grid level.

Parameters

user	The UserCtx for which the DMs will be created. Its IM, JM, KM fields must be pre-populated.
coarse_user	The UserCtx of the next-coarser grid level, or NULL if user is the coarsest level.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 123 of file grid.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
 
    DMBoundaryType xperiod = (simCtx->i_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    DMBoundaryType yperiod = (simCtx->j_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    DMBoundaryType zperiod = (simCtx->k_periodic) ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE;
    PetscInt stencil_width = (simCtx->i_periodic || simCtx->j_periodic || simCtx->k_periodic) ? 3:2; // Stencil width is 2 in the legacy code
 
    PetscInt *lx = NULL, *ly = NULL, *lz = NULL;
    PetscInt m, n, p;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (coarse_user) {
        // --- This is a FINE grid; it must be aligned with the COARSE grid ---
        LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: [Aligning DM] for block %d level %d (size %dx%dx%d) with level %d\n", simCtx->rank, user->_this, user->thislevel, user->IM, user->JM, user->KM, coarse_user->thislevel);
 
        DMDAGetInfo(coarse_user->da, NULL, NULL, NULL, NULL, &m, &n, &p, NULL, NULL, NULL, NULL, NULL, NULL);
        LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d:   Coarse grid processor decomposition is %d x %d x %d\n", simCtx->rank, m, n, p);
 
        // This is the core logic from MGDACreate to ensure processor alignment.
        PetscInt *lx_contrib, *ly_contrib, *lz_contrib;
        ierr = PetscMalloc3(m, &lx_contrib, n, &ly_contrib, p, &lz_contrib); CHKERRQ(ierr);
        ierr = PetscMemzero(lx_contrib, m * sizeof(PetscInt)); CHKERRQ(ierr);
        ierr = PetscMemzero(ly_contrib, n * sizeof(PetscInt)); CHKERRQ(ierr);
        ierr = PetscMemzero(lz_contrib, p * sizeof(PetscInt)); CHKERRQ(ierr);
 
        DMDALocalInfo info;
        DMDAGetLocalInfo(coarse_user->da, &info);
        PetscInt xs = info.xs, xe = info.xs + info.xm, mx = info.mx;
        PetscInt ys = info.ys, ye = info.ys + info.ym, my = info.my;
        PetscInt zs = info.zs, ze = info.zs + info.zm, mz = info.mz;
        
        PetscMPIInt rank;
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
        PetscInt proc_i = rank % m;
        PetscInt proc_j = (rank / m) % n;
        PetscInt proc_k = rank / (m * n);
 
        // --- X-Direction Logic (Identical to MGDACreate) ---
        if (user->isc) lx_contrib[proc_i] = (xe - xs);
        else {
            if (m == 1) lx_contrib[0] = user->IM + 1;
            else if (xs == 0) lx_contrib[0] = 2 * xe - 1;
            else if (xe == mx) lx_contrib[proc_i] = user->IM + 1 - (2 * xs - 1);
            else lx_contrib[proc_i] = (xe - xs) * 2;
        }
 
        // --- Y-Direction Logic (Identical to MGDACreate) ---
        if (user->jsc) ly_contrib[proc_j] = (ye - ys);
        else {
            if (n == 1) ly_contrib[0] = user->JM + 1;
            else if (ys == 0) ly_contrib[0] = 2 * ye - 1;
            else if (ye == my) ly_contrib[proc_j] = user->JM + 1 - (2 * ys - 1);
            else ly_contrib[proc_j] = (ye - ys) * 2;
        }
 
        // --- Z-Direction Logic (Identical to MGDACreate) ---
        if (user->ksc) lz_contrib[proc_k] = (ze - zs);
        else {
            if (p == 1) lz_contrib[0] = user->KM + 1;
            else if (zs == 0) lz_contrib[0] = 2 * ze - 1;
            else if (ze == mz) lz_contrib[proc_k] = user->KM + 1 - (2 * zs - 1);
            else lz_contrib[proc_k] = (ze - zs) * 2;
        }
        LOG_ALLOW_SYNC(LOCAL, LOG_VERBOSE, "Rank %d:   Calculated this rank's node contribution to fine grid: lx=%d, ly=%d, lz=%d\n", simCtx->rank, lx_contrib[proc_i], ly_contrib[proc_j], lz_contrib[proc_k]);
 
        // Allocate the final distribution arrays and Allreduce to get the global distribution
        ierr = PetscMalloc3(m, &lx, n, &ly, p, &lz); CHKERRQ(ierr);
        ierr = MPI_Allreduce(lx_contrib, lx, m, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(ly_contrib, ly, n, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(lz_contrib, lz, p, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        ierr = PetscFree3(lx_contrib, ly_contrib, lz_contrib); CHKERRQ(ierr);
 
    } else {
        // --- CASE 2: This is the COARSEST grid; use default or user-specified decomposition ---
        if(simCtx->exec_mode == EXEC_MODE_SOLVER){
 
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating coarsest DM for block %d level %d (size %dx%dx%d)\n", simCtx->rank, user->_this, user->thislevel, user->IM, user->JM, user->KM);
            m = simCtx->da_procs_x;
            n = simCtx->da_procs_y;
            p = simCtx->da_procs_z;
        
        } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
 
          LOG_ALLOW(GLOBAL,LOG_ERROR,"Currently Only Single Rank is supported. \n");
          
          m = n = p = PETSC_DECIDE;
            
        } 
        // lx, ly, lz are NULL, so DMDACreate3d will use the m,n,p values.
    }
 
    // --- Create the DMDA for the current UserCtx ---
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d:   Calling DMDACreate3d...\n", simCtx->rank);
    ierr = DMDACreate3d(PETSC_COMM_WORLD, xperiod, yperiod, zperiod, DMDA_STENCIL_BOX,
                          user->IM + 1, user->JM + 1, user->KM + 1,
                          m, n, p,
                          1, stencil_width, lx, ly, lz, &user->da); CHKERRQ(ierr);
    
    if (coarse_user) {
        ierr = PetscFree3(lx, ly, lz); CHKERRQ(ierr);
    }
    
    // --- Standard DM setup applicable to all levels ---
    ierr = DMSetUp(user->da); CHKERRQ(ierr);
    ierr = DMGetCoordinateDM(user->da, &user->fda); CHKERRQ(ierr);
    ierr = DMDASetUniformCoordinates(user->da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(user->da, &user->info); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d:   DM creation for block %d level %d complete.\n", simCtx->rank, user->_this, user->thislevel);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

InitializeAllGridDMs()

PetscErrorCode InitializeAllGridDMs

(

SimCtx *

simCtx

)

Orchestrates the creation of DMDA objects for every block and multigrid level.

This function systematically builds the entire DMDA hierarchy. It first calculates the dimensions (IM, JM, KM) for all coarse grids based on the finest grid's dimensions and the semi-coarsening flags. It then iterates from the coarsest to the finest level, calling a powerful helper function (InitializeSingleGridDM) to create the DMs for each block, ensuring that finer grids are properly aligned with their coarser parents for multigrid efficiency.

Parameters

simCtx

The master SimCtx, containing the configured UserCtx hierarchy.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 260 of file grid.c.

{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx = usermg->mgctx;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "Pre-scanning BCs to identify domain periodicity.\n");
    ierr = DeterminePeriodicity(simCtx); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Creating DMDA objects for all levels and blocks...\n");
 
    // --- Part 1: Calculate Coarse Grid Dimensions & VALIDATE ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating and validating coarse grid dimensions...\n");
    for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *user_coarse = &mgctx[level].user[bi];
            UserCtx *user_fine   = &mgctx[level + 1].user[bi];
 
            user_coarse->IM = user_fine->isc ? user_fine->IM : (user_fine->IM + 1) / 2;
            user_coarse->JM = user_fine->jsc ? user_fine->JM : (user_fine->JM + 1) / 2;
            user_coarse->KM = user_fine->ksc ? user_fine->KM : (user_fine->KM + 1) / 2;
 
            LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Block %d, Level %d dims calculated: %d x %d x %d\n",
                      simCtx->rank, bi, level, user_coarse->IM, user_coarse->JM, user_coarse->KM);
 
            // Validation check from legacy MGDACreate to ensure coarsening is possible
            PetscInt check_i = user_coarse->IM * (2 - user_coarse->isc) - (user_fine->IM + 1 - user_coarse->isc);
            PetscInt check_j = user_coarse->JM * (2 - user_coarse->jsc) - (user_fine->JM + 1 - user_coarse->jsc);
            PetscInt check_k = user_coarse->KM * (2 - user_coarse->ksc) - (user_fine->KM + 1 - user_coarse->ksc);
 
            if (check_i + check_j + check_k != 0) {
          //      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
          //           "Grid at level %d, block %d cannot be coarsened from %dx%dx%d to %dx%dx%d with the given semi-coarsening flags. Check grid dimensions.",
              //          level, bi, user_fine->IM, user_fine->JM, user_fine->KM, user_coarse->IM, user_coarse->JM, user_coarse->KM);
          LOG(GLOBAL,LOG_WARNING,"WARNING: Grid at level %d, block %d can't be consistently coarsened further.\n", level, bi);
            }
        }
    }
 
    // --- Part 2: Create DMs from Coarse to Fine for each Block ---
    for (PetscInt bi = 0; bi < nblk; bi++) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "--- Creating DMs for Block %d ---\n", bi);
 
        // Create the coarsest level DM first (passing NULL for the coarse_user)
        ierr = InitializeSingleGridDM(&mgctx[0].user[bi], NULL); CHKERRQ(ierr);
 
        // Create finer level DMs, passing the next-coarser context for alignment
        for (PetscInt level = 1; level < usermg->mglevels; level++) {
            ierr = InitializeSingleGridDM(&mgctx[level].user[bi], &mgctx[level-1].user[bi]); CHKERRQ(ierr);
        }
    }
 
    // --- Optional: View the finest DM for debugging verification ---
    if (get_log_level() >= LOG_DEBUG) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "--- Viewing Finest DMDA (Level %d, Block 0) ---\n", usermg->mglevels - 1);
        ierr = DMView(mgctx[usermg->mglevels - 1].user[0].da, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "DMDA object creation complete.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetFinestLevelCoordinates()

static PetscErrorCode SetFinestLevelCoordinates ( UserCtx *  user )

static

A router that populates the coordinates for a single finest-level DMDA.

This function orchestrates the coordinate setting for one block. It checks the global generate_grid flag and calls the appropriate helper for either programmatic generation or reading from a file.

After the local coordinate vector is populated by a helper, this function performs the necessary DMLocalToGlobal and DMGlobalToLocal scatters to ensure that the ghost node coordinate values are correctly communicated and updated across all MPI ranks.

Parameters

user	The UserCtx for a specific block on the finest level.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 420 of file grid.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting finest level coordinates for block %d...\n", simCtx->rank, user->_this);
 
    if (simCtx->generate_grid) {
        ierr = GenerateAndSetCoordinates(user); CHKERRQ(ierr);
    } else {
 
      FILE *grid_file_handle = NULL;
      // Only Rank 0 opens the file.
      if (simCtx->rank == 0) {
    grid_file_handle = fopen(simCtx->grid_file, "r");
    if (!grid_file_handle) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file: %s", simCtx->grid_file);
        
    // Now, on Rank 0, we skip the entire header section once.
    // This is the logic from your modern code's AssignGridCoordinates.
    PetscInt headerLines = simCtx->block_number + 2; // 1 for nblk, plus one for each block's dims
    char dummy_buffer[2048];
    for (PetscInt s = 0; s < headerLines; ++s) {
      if (!fgets(dummy_buffer, sizeof(dummy_buffer), grid_file_handle)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unexpected EOF while skipping grid header");
      }
    }
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank 0: Skipped %d header lines, now at coordinate data.\n", headerLines);
      }
 
      // We now call the coordinate reader, passing the file handle.
      // It's responsible for reading its block's data and broadcasting.
      ierr = ReadAndSetCoordinates(user, grid_file_handle); CHKERRQ(ierr);
        
      // Only Rank 0, which opened the file, should close it.
      if (simCtx->rank == 0) {
    fclose(grid_file_handle);
      }
    }
    
    // After populating the local coordinate vector, we must perform a
    // Local-to-Global and then Global-to-Local scatter to correctly
    // populate the ghost node coordinates across process boundaries.
    Vec gCoor, lCoor;
    ierr = DMGetCoordinates(user->da, &gCoor); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(user->da, &lCoor); CHKERRQ(ierr);
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Scattering coordinates to update ghost nodes for block %d...\n", simCtx->rank, user->_this);
    ierr = DMLocalToGlobalBegin(user->fda, lCoor, INSERT_VALUES, gCoor); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(user->fda, lCoor, INSERT_VALUES, gCoor); CHKERRQ(ierr);
    
    ierr = DMGlobalToLocalBegin(user->fda, gCoor, INSERT_VALUES, lCoor); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(user->fda, gCoor, INSERT_VALUES, lCoor); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GenerateAndSetCoordinates()

static PetscErrorCode GenerateAndSetCoordinates ( UserCtx *  user )

static

Programmatically generates and sets grid coordinates based on user parameters.

This function populates the local coordinate vector of the provided UserCtx using the geometric properties (Min_X, Max_X, rx, etc.) that were parsed from command-line options. It supports non-linear grid stretching.

Parameters

user	The UserCtx for a specific block.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

DEBUG: This verifies the presence of a last "unphysical" layer of coordinates.

Definition at line 533 of file grid.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Cmpnts       ***coor;
    Vec            lCoor;
    
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Generating coordinates for block %d...\n", user->simCtx->rank, user->_this);
    
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(user->da, &lCoor); CHKERRQ(ierr);
 
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Local Info for block %d - X range - [%d,%d], Y range - [%d,%d], Z range - [%d,%d]\n",
              user->simCtx->rank, user->_this, xs, xe, ys, ye, zs, ze);
    
    LOG_ALLOW_SYNC(LOCAL, LOG_TRACE, "Rank %d: Local Info for block %d - X domain - [%.4f,%.4f], Y range - [%.4f,%.4f], Z range - [%.4f,%.4f]\n",
              user->simCtx->rank, user->_this, user->Min_X,user->Max_X,user->Min_Y,user->Max_Y,user->Min_Z,user->Max_Z);          
              
    ierr = VecSet(lCoor, 0.0); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, lCoor, &coor); CHKERRQ(ierr);
    
    PetscReal Lx = user->Max_X - user->Min_X;
    PetscReal Ly = user->Max_Y - user->Min_Y;
    PetscReal Lz = user->Max_Z - user->Min_Z;
 
    // Loop over the local nodes, including ghost nodes, owned by this process.
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                if(k < user->KM && j < user->JM && i < user->IM){
                    coor[k][j][i].x = user->Min_X + ComputeStretchedCoord(i, user->IM, Lx, user->rx);
                    coor[k][j][i].y = user->Min_Y + ComputeStretchedCoord(j, user->JM, Ly, user->ry);
                    coor[k][j][i].z = user->Min_Z + ComputeStretchedCoord(k, user->KM, Lz, user->rz);
                }
            }  
        }
    }
 
    /// DEBUG: This verifies the presence of a last "unphysical" layer of coordinates.
    /*
    PetscInt KM = user->KM;
    for (PetscInt j = ys; j < ye; j++){
        for(PetscInt i = xs; i < xe; i++){
            LOG_ALLOW(GLOBAL,LOG_DEBUG,"coor[%d][%d][%d].(x,y,z) = %le,%le,%le",KM,j,i,coor[KM][j][i].x,coor[KM][j][i].y,coor[KM][j][i].z);
        }
    }
    */
 
 
 
    ierr = DMDAVecRestoreArray(user->fda, lCoor, &coor); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadAndSetCoordinates()

static PetscErrorCode ReadAndSetCoordinates ( UserCtx *  user, FILE *  fd  )

static

Reads physical coordinates from a file and populates the DMDA for a specific block.

This function handles the collective read of an interleaved (X Y Z per line) multi-block grid file. It assumes the file header (nblk, dimensions) has already been processed by ReadGridFile.

The process is robust for parallel execution:

1. Rank 0 opens the grid file.
2. It intelligently skips past the header section and the coordinate data for all blocks preceding the current block being processed (user->_this).
3. It then reads the entire coordinate data for the current block into a single contiguous buffer gc.
4. This global buffer gc is broadcast to all other MPI ranks.
5. Each rank then loops through its local (owned + ghost) node indices, calculates the corresponding index in the global gc array, and copies the (x,y,z) coordinates into its local PETSc coordinate vector.

Parameters

user	The UserCtx for a specific block. Its _this field must be set, and its IM, JM, KM fields must be correctly pre-populated.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 621 of file grid.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
    PetscMPIInt    rank = simCtx->rank;
    PetscInt       block_index = user->_this;
    PetscInt       IM = user->IM, JM = user->JM, KM = user->KM;
    DMDALocalInfo  info;
    Cmpnts       ***coor;
    Vec            lCoor;
    PetscReal      *gc = NULL; // Global coordinate buffer, allocated on all ranks
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Reading interleaved coordinates from file for block %d...\n",
              simCtx->rank, block_index);
 
    // 1. Allocate the buffer on ALL ranks to receive the broadcast data.
    // PetscInt n_nodes = (IM + 1) * (JM + 1) * (KM + 1);
    PetscInt n_nodes = (IM) * (JM) * (KM);
    ierr = PetscMalloc1(3 * n_nodes, &gc); CHKERRQ(ierr);
 
    // 2. Only Rank 0 opens the file and reads the data.
    if (rank == 0) {
        if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Recieved a NULL file handle.\n");
                
        // Read the coordinate data for the CURRENT block.
        for (PetscInt k = 0; k < KM; k++) {
            for (PetscInt j = 0; j < JM; j++) {
                for (PetscInt i = 0; i < IM; i++) {
                    PetscInt base_index = 3 * ((k * (JM) + j) * (IM) + i);
                    if (fscanf(fd, "%le %le %le\n", &gc[base_index], &gc[base_index + 1], &gc[base_index + 2]) != 3) {
                        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Error reading coordinates for node (i,j,k)=(%d,%d,%d) in block %d", i, j, k, block_index);
                    }
                }
            }
        }
        
    }
 
    // 3. Broadcast the coordinate block for the current block to all other processes.
    ierr = MPI_Bcast(gc, 3 * n_nodes, MPIU_REAL, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // 4. Each rank populates its local portion of the coordinate vector.
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(user->da, &lCoor); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, lCoor, &coor); CHKERRQ(ierr);
 
    for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
        for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
            for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
          if(k< KM && j < JM && i < IM){
                PetscInt base_idx = 3 * ((k * (JM) + j) * (IM) + i);
                coor[k][j][i].x = gc[base_idx];
                coor[k][j][i].y = gc[base_idx + 1];
                coor[k][j][i].z = gc[base_idx + 2];
          }
        }
        }
    }
    
    // 5. Clean up and restore.
    ierr = DMDAVecRestoreArray(user->fda, lCoor, &coor); CHKERRQ(ierr);
    ierr = PetscFree(gc); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

RestrictCoordinates()

static PetscErrorCode RestrictCoordinates ( UserCtx *  coarse_user, UserCtx *  fine_user  )

static

Populates coarse grid coordinates by restricting from a fine grid.

This function is a direct adaptation of the coordinate restriction loop from the legacy MG_Initial function. It ensures that the physical location of a coarse grid node is identical to its corresponding parent node on the fine grid. The mapping from coarse to fine index (i -> ih) is determined by the semi-coarsening flags (isc, jsc, ksc) stored in the UserCtx.

Parameters

coarse_user	The UserCtx for the coarse grid (destination).
fine_user	The UserCtx for the fine grid (source).

Returns

PetscErrorCode

Definition at line 708 of file grid.c.

{
    PetscErrorCode ierr;
    Vec            c_lCoor, f_lCoor; // Coarse and Fine local coordinate vectors
    Cmpnts       ***c_coor;
    const Cmpnts ***f_coor; // Use const for read-only access
    DMDALocalInfo  c_info;
    PetscInt       ih, jh, kh; // Fine-grid indices corresponding to coarse-grid i,j,k
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Restricting coords from level %d to level %d for block %d\n",
                   fine_user->simCtx->rank, fine_user->thislevel, coarse_user->thislevel, coarse_user->_this);
 
    ierr = DMDAGetLocalInfo(coarse_user->da, &c_info); CHKERRQ(ierr);
 
    ierr = DMGetCoordinatesLocal(coarse_user->da, &c_lCoor); CHKERRQ(ierr);
    ierr = DMGetCoordinatesLocal(fine_user->da, &f_lCoor); CHKERRQ(ierr);
    
    ierr = DMDAVecGetArray(coarse_user->fda, c_lCoor, &c_coor); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(fine_user->fda, f_lCoor, &f_coor); CHKERRQ(ierr);
 
    // Get the local owned range of the coarse grid.
    PetscInt xs = c_info.xs, xe = c_info.xs + c_info.xm;
    PetscInt ys = c_info.ys, ye = c_info.ys + c_info.ym;
    PetscInt zs = c_info.zs, ze = c_info.zs + c_info.zm;
 
    // Get the global dimensions of the coarse grid.
    PetscInt mx = c_info.mx, my = c_info.my, mz = c_info.mz;
 
    // If this process owns the maximum boundary node, contract the loop by one
    // to prevent the index doubling `2*i` from going out of bounds.
    // This is also ensuring we do not manipulate the unphysical layer of coors present in the finest level.
    if (xe == mx) xe--;
    if (ye == my) ye--;
    if (ze == mz) ze--;
    
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                // Determine the corresponding parent node index on the FINE grid,
                // respecting the semi-coarsening flags of the FINE grid's UserCtx.
                ih = coarse_user->isc ? i : 2 * i;
                jh = coarse_user->jsc ? j : 2 * j;
        kh = coarse_user->ksc ? k : 2 * k;
 
        //  LOG_ALLOW(GLOBAL,LOG_DEBUG," [kh][ih][jh] = %d,%d,%d  - k,j,i = %d,%d,%d.\n",kh,jh,ih,k,j,i);
 
                c_coor[k][j][i] = f_coor[kh][jh][ih];
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(coarse_user->fda, c_lCoor, &c_coor); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(fine_user->fda, f_lCoor, &f_coor); CHKERRQ(ierr);
 
    // After populating the local portion, scatter to update ghost regions.
    Vec c_gCoor;
    ierr = DMGetCoordinates(coarse_user->da, &c_gCoor); CHKERRQ(ierr);
    ierr = DMLocalToGlobalBegin(coarse_user->fda, c_lCoor, INSERT_VALUES, c_gCoor); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(coarse_user->fda, c_lCoor, INSERT_VALUES, c_gCoor); CHKERRQ(ierr);
    ierr = DMGlobalToLocalBegin(coarse_user->fda, c_gCoor, INSERT_VALUES, c_lCoor); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(coarse_user->fda, c_gCoor, INSERT_VALUES, c_lCoor); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

AssignAllGridCoordinates()

PetscErrorCode AssignAllGridCoordinates

(

SimCtx *

simCtx

)

Orchestrates the assignment of physical coordinates to all DMDA objects.

This function manages the entire process of populating the coordinate vectors for every DMDA across all multigrid levels and blocks. It follows a two-part strategy that is essential for multigrid methods:

1. Populate Finest Level: It first loops through each block and calls a helper (SetFinestLevelCoordinates) to set the physical coordinates for the highest-resolution grid (the finest multigrid level).
2. Restrict to Coarser Levels: It then iterates downwards from the finest level, calling a helper (RestrictCoordinates) to copy the coordinate values from the fine grid nodes to their corresponding parent nodes on the coarser grids. This ensures all levels represent the exact same geometry.

Parameters

simCtx

The master SimCtx, containing the configured UserCtx hierarchy.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 356 of file grid.c.

{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Assigning physical coordinates to all grid DMs...\n");
 
    // --- Part 1: Populate the Finest Grid Level ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Setting coordinates for the finest grid level (%d)...\n", usermg->mglevels - 1);
    for (PetscInt bi = 0; bi < nblk; bi++) {
        UserCtx *fine_user = &usermg->mgctx[usermg->mglevels - 1].user[bi];
        ierr = SetFinestLevelCoordinates(fine_user); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL,LOG_TRACE,"The Finest level coordinates for block %d have been set.\n",bi);
        if(get_log_level()==LOG_VERBOSE && is_function_allowed(__FUNCT__)==true){
            ierr = LOG_FIELD_MIN_MAX(fine_user,"Coordinates");
        }
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finest level coordinates have been set for all blocks.\n");
 
    // --- Part 2: Restrict Coordinates to Coarser Levels ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Restricting coordinates to coarser grid levels...\n");
    for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *coarse_user = &usermg->mgctx[level].user[bi];
            UserCtx *fine_user   = &usermg->mgctx[level + 1].user[bi];
            ierr = RestrictCoordinates(coarse_user, fine_user); CHKERRQ(ierr);
 
            LOG_ALLOW(GLOBAL,LOG_TRACE,"Coordinates restricted to block %d level %d.\n",bi,level);
            if(get_log_level==LOG_VERBOSE && is_function_allowed(__FUNCT__)==true) {
                ierr = LOG_FIELD_MIN_MAX(coarse_user,"Coordinates");
            }
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Physical coordinates assigned to all grid levels and blocks.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeStretchedCoord()

static PetscReal ComputeStretchedCoord ( PetscInt  i, PetscInt  N, PetscReal  L, PetscReal  r  )

inlinestatic

Computes a stretched coordinate along one dimension.

This function computes a coordinate based on a geometric stretching ratio. If the ratio (r) is 1.0, a uniform distribution is used: x(i) = L * (i/N)

If r != 1.0, a geometric stretching is applied: x(i) = L * [ (r^(i/N) - 1 ) / (r - 1) ]

Here:

• i : The current index along the dimension.
• N : The total number of divisions along that dimension.
• L : The length of the domain along that dimension.
• r : The stretching ratio. r > 1.0 stretches the grid in a geometric fashion increasing spacing away from the start, whereas 0 < r < 1.0 would cluster points near the start.

Parameters

[in]	i	Current index (0 <= i <= N).
[in]	N	Number of segments along the dimension.
[in]	L	Total length of the domain.
[in]	r	Stretching ratio.

Returns

PetscReal The computed coordinate at index i.

Note

This function does not return a PetscErrorCode because it does not allocate memory or call PETSc routines that can fail. It is just a helper calculation function.

Definition at line 510 of file grid.c.

{
    if (N <=1) return 0.0;
    PetscReal fraction = (PetscReal)i / ((PetscReal)N - 1.0);
    if (PetscAbsReal(r - 1.0) < 1.0e-9) { // Use a tolerance for float comparison
        return L * fraction;
    } else {
        return L * (PetscPowReal(r, fraction) - 1.0) / (r - 1.0);
    }
}

Here is the caller graph for this function:

ComputeLocalBoundingBox()

PetscErrorCode ComputeLocalBoundingBox	(	UserCtx *	user,
		BoundingBox *	localBBox
	)

Computes the local bounding box of the grid on the current process.

This function calculates the minimum and maximum coordinates (x, y, z) of the local grid points owned by the current MPI process. It iterates over the local portion of the grid, examines each grid point's coordinates, and updates the minimum and maximum values accordingly.

The computed bounding box is stored in the provided localBBox structure, and the user->bbox field is also updated with this bounding box for consistency within the user context.

Parameters

[in]	user	Pointer to the user-defined context containing grid information. This context must be properly initialized before calling this function.
[out]	localBBox	Pointer to the BoundingBox structure where the computed local bounding box will be stored. The structure should be allocated by the caller.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Definition at line 800 of file grid.c.

{
    PetscErrorCode ierr;
    PetscInt i, j, k;
    PetscMPIInt rank;
    PetscInt xs, ys, zs, xe, ye, ze;
    DMDALocalInfo info;
    Vec coordinates;
    Cmpnts ***coordArray;
    Cmpnts minCoords, maxCoords;
 
    PetscFunctionBeginUser;
    
    PROFILE_FUNCTION_BEGIN;
 
    // Start of function execution
    LOG_ALLOW(GLOBAL, LOG_INFO, "Entering the function.\n");
 
    // Validate input Pointers
    if (!user) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Input 'user' Pointer is NULL.\n");
        return PETSC_ERR_ARG_NULL;
    }
    if (!localBBox) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Output 'localBBox' Pointer is NULL.\n");
        return PETSC_ERR_ARG_NULL;
    }
 
    // Get MPI rank
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // Get the local coordinates vector from the DMDA
    ierr = DMGetCoordinatesLocal(user->da, &coordinates);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error getting local coordinates vector.\n");
        return ierr;
    }
 
    if (!coordinates) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Coordinates vector is NULL.\n");
        return PETSC_ERR_ARG_NULL;
    }
 
    // Access the coordinate array for reading
    ierr = DMDAVecGetArrayRead(user->fda, coordinates, &coordArray);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error accessing coordinate array.\n");
        return ierr;
    }
 
    // Get the local grid information (indices and sizes)
    ierr = DMDAGetLocalInfo(user->da, &info);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error getting DMDA local info.\n");
        return ierr;
    }
 
    
    xs = info.gxs; xe = xs + info.gxm;
    ys = info.gys; ye = ys + info.gym;
    zs = info.gzs; ze = zs + info.gzm;
    
    /*
    xs = info.xs; xe = xs + info.xm;
    ys = info.ys; ye = ys + info.ym;
    zs = info.zs; ze = zs + info.zm;
    */
    
    // Initialize min and max coordinates with extreme values
    minCoords.x = minCoords.y = minCoords.z = PETSC_MAX_REAL;
    maxCoords.x = maxCoords.y = maxCoords.z = PETSC_MIN_REAL;
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Grid indices (Including Ghosts): xs=%d, xe=%d, ys=%d, ye=%d, zs=%d, ze=%d.\n",rank, xs, xe, ys, ye, zs, ze);
 
    // Iterate over the local grid to find min and max coordinates
    for (k = zs; k < ze; k++) {
        for (j = ys; j < ye; j++) {
            for (i = xs; i < xe; i++) {
                // Only consider nodes within the physical domain.
                if(i < user->IM && j < user->JM && k < user->KM){
                    Cmpnts coord = coordArray[k][j][i];
 
                    // Update min and max coordinates
                    if (coord.x < minCoords.x) minCoords.x = coord.x;
                    if (coord.y < minCoords.y) minCoords.y = coord.y;
                    if (coord.z < minCoords.z) minCoords.z = coord.z;
 
                    if (coord.x > maxCoords.x) maxCoords.x = coord.x;
                    if (coord.y > maxCoords.y) maxCoords.y = coord.y;
                    if (coord.z > maxCoords.z) maxCoords.z = coord.z;
                }    
            }
        }
    }
 
 
    // Add tolerance to bboxes.
    minCoords.x =  minCoords.x - BBOX_TOLERANCE;
    minCoords.y =  minCoords.y - BBOX_TOLERANCE;
    minCoords.z =  minCoords.z - BBOX_TOLERANCE;
 
    maxCoords.x =  maxCoords.x + BBOX_TOLERANCE;
    maxCoords.y =  maxCoords.y + BBOX_TOLERANCE;
    maxCoords.z =  maxCoords.z + BBOX_TOLERANCE;
 
    LOG_ALLOW(LOCAL,LOG_DEBUG," Tolerance added to the limits: %.8e .\n",(PetscReal)BBOX_TOLERANCE);
       
    // Log the computed min and max coordinates
     LOG_ALLOW(LOCAL, LOG_INFO,"[Rank %d] Bounding Box Ranges = X[%.6f, %.6f], Y[%.6f,%.6f], Z[%.6f, %.6f].\n",rank,minCoords.x, maxCoords.x,minCoords.y, maxCoords.y, minCoords.z, maxCoords.z);
 
 
    
    // Restore the coordinate array
    ierr = DMDAVecRestoreArrayRead(user->fda, coordinates, &coordArray);
    if (ierr) {
        LOG_ALLOW(LOCAL, LOG_ERROR, "Error restoring coordinate array.\n");
        return ierr;
    }
 
    // Set the local bounding box
    localBBox->min_coords = minCoords;
    localBBox->max_coords = maxCoords;
 
    // Update the bounding box inside the UserCtx for consistency
    user->bbox = *localBBox;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Exiting the function successfully.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

GatherAllBoundingBoxes()

PetscErrorCode GatherAllBoundingBoxes	(	UserCtx *	user,
		BoundingBox **	allBBoxes
	)

Gathers local bounding boxes from all MPI processes to rank 0.

Each rank computes its local bounding box, then all ranks participate in an MPI_Gather to send their BoundingBox to rank 0. Rank 0 allocates the result array and returns it via allBBoxes.

Parameters

[in]	user	Pointer to UserCtx (must be non-NULL).
[out]	allBBoxes	On rank 0, receives malloc’d array of size size. On other ranks, set to NULL.

Returns

PetscErrorCode

Definition at line 948 of file grid.c.

{
  PetscErrorCode ierr;
  PetscMPIInt    rank, size;
  BoundingBox    *bboxArray = NULL;
  BoundingBox     localBBox;
 
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  /* Validate */
  if (!user || !allBBoxes) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
                                    "GatherAllBoundingBoxes: NULL pointer");
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRMPI(ierr);
 
  /* Compute local bbox */
  ierr = ComputeLocalBoundingBox(user, &localBBox); CHKERRQ(ierr);
 
  /* Ensure everyone is synchronized before the gather */
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: about to MPI_Gather(localBBox)\n", rank);
 
  /* Allocate on root */
  if (rank == 0) {
    bboxArray = (BoundingBox*)malloc(size * sizeof(BoundingBox));
    if (!bboxArray) SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_MEM,
                                "GatherAllBoundingBoxes: malloc failed");
  }
 
  /* Collective: every rank must call */
  ierr = MPI_Gather(&localBBox, sizeof(BoundingBox), MPI_BYTE,
                    bboxArray, sizeof(BoundingBox), MPI_BYTE,
                    0, PETSC_COMM_WORLD);
  CHKERRMPI(ierr);
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: completed MPI_Gather(localBBox)\n", rank);
 
  /* Return result */
  if (rank == 0) {
    *allBBoxes = bboxArray;
  } else {
    *allBBoxes = NULL;
  }
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BroadcastAllBoundingBoxes()

PetscErrorCode BroadcastAllBoundingBoxes	(	UserCtx *	user,
		BoundingBox **	bboxlist
	)

Broadcasts the bounding box list from rank 0 to all ranks.

Broadcasts the bounding box information collected on rank 0 to all other ranks.

After GatherAllBoundingBoxes, rank 0 has an array of size boxes. This routine makes sure every rank ends up with its own malloc’d copy.

Parameters

[in]	user	Pointer to UserCtx (unused here, but kept for signature).
[in,out]	bboxlist	On entry: rank 0’s array; on exit: every rank’s array.

Returns

PetscErrorCode

Definition at line 1016 of file grid.c.

{
  PetscErrorCode ierr;
  PetscMPIInt    rank, size;
 
  PetscFunctionBeginUser;
 
  PROFILE_FUNCTION_BEGIN;
 
  if (!bboxlist) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
                          "BroadcastAllBoundingBoxes: NULL pointer");
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRMPI(ierr);
 
  /* Non-root ranks must allocate before the Bcast */
  if (rank != 0) {
    *bboxlist = (BoundingBox*)malloc(size * sizeof(BoundingBox));
    if (!*bboxlist) SETERRABORT(PETSC_COMM_WORLD, PETSC_ERR_MEM,
                                "BroadcastAllBoundingBoxes: malloc failed");
  }
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: about to MPI_Bcast(%d boxes)\n", rank, size);
 
  /* Collective: every rank must call */
  ierr = MPI_Bcast(*bboxlist, size * sizeof(BoundingBox), MPI_BYTE,
                   0, PETSC_COMM_WORLD);
  CHKERRMPI(ierr);
 
  MPI_Barrier(PETSC_COMM_WORLD);
  LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE,
    "Rank %d: completed MPI_Bcast(%d boxes)\n", rank, size);
 
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the caller graph for this function:

CalculateInletProperties()

PetscErrorCode CalculateInletProperties

(

UserCtx *

user

)

Calculates the center and area of the primary INLET face.

This function identifies the primary INLET face from the boundary face configurations, computes its geometric center and total area using a generic utility function, and stores these results in the simulation context.

Parameters

user	Pointer to the UserCtx containing boundary face information.

Returns

PetscErrorCode

Definition at line 1069 of file grid.c.

{
    PetscErrorCode ierr;
    BCFace         inlet_face_id = -1;
    PetscBool      inlet_found = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // 1. Identify the primary inlet face from the configuration
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == INLET) {
            inlet_face_id = user->boundary_faces[i].face_id;
            inlet_found = PETSC_TRUE;
            break; // Use the first inlet found
        }
    }
 
    if (!inlet_found) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "No INLET face found. Skipping inlet center calculation.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    
    Cmpnts inlet_center;
    PetscReal inlet_area;
 
    // 2. Call the generic utility to compute  the center and area of any face.
    ierr = CalculateFaceCenterAndArea(user,inlet_face_id,&inlet_center,&inlet_area); CHKERRQ(ierr);
 
    // 3. Store results in the SimCtx
    user->simCtx->CMx_c = inlet_center.x;
    user->simCtx->CMy_c = inlet_center.y;
    user->simCtx->CMz_c = inlet_center.z;
    user->simCtx->AreaInSum = inlet_area;
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Rank[%d] Inlet Center: (%.6f, %.6f, %.6f), Area: %.6f\n",
              inlet_center.x, inlet_center.y, inlet_center.z, inlet_area);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
 
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateOutletProperties()

PetscErrorCode CalculateOutletProperties

(

UserCtx *

user

)

Calculates the center and area of the primary OUTLET face.

This function identifies the primary OUTLET face from the boundary face configurations, computes its geometric center and total area using a generic utility function, and stores these results in the simulation context.

Parameters

user	Pointer to the UserCtx containing boundary face information.

Returns

PetscErrorCode

Definition at line 1125 of file grid.c.

{
    PetscErrorCode ierr;
    BCFace         outlet_face_id = -1;
    PetscBool      outlet_found = PETSC_FALSE;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    // 1. Identify the primary outlet face from the configuration
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == OUTLET) {
            outlet_face_id = user->boundary_faces[i].face_id;
            outlet_found = PETSC_TRUE;
            break; // Use the first outlet found
        }
    }
    if (!outlet_found) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "No OUTLET face found. Skipping outlet center calculation.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    PetscReal outlet_area;
    Cmpnts outlet_center;
    // 2. Call the generic utility to compute  the center and area of any face
    ierr = CalculateFaceCenterAndArea(user,outlet_face_id,&outlet_center,&outlet_area); CHKERRQ(ierr);
    // 3. Store results in the SimCtx
    user->simCtx->AreaOutSum = outlet_area;
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Outlet Center: (%.6f, %.6f, %.6f), Area: %.6f\n",
              outlet_center.x, outlet_center.y, outlet_center.z, outlet_area);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateFaceCenterAndArea()

PetscErrorCode CalculateFaceCenterAndArea	(	UserCtx *	user,
		BCFace	face_id,
		Cmpnts *	face_center,
		PetscReal *	face_area
	)

Calculates the geometric center and total area of a specified boundary face.

This function computes two key properties of a boundary face in the computational domain:

1. Geometric Center: The average (x,y,z) position of all physical nodes on the face
2. Total Area: The sum of face area vector magnitudes from all non-solid cells adjacent to the face

Indexing Architecture

The solver uses different indexing conventions for different field types:

Node-Centered Fields (Coordinates):

• Direct indexing: Node n stored at coor[n]
• For mx=26: Physical nodes [0-24], Dummy at [25]
• For mz=98: Physical nodes [0-96], Dummy at [97]

Face-Centered Fields (Metrics: csi, eta, zet):

• Direct indexing: Face n stored at csi/eta/zet[n]
• For mx=26: Physical faces [0-24], Dummy at [25]
• For mz=98: Physical faces [0-96], Dummy at [97]
• Face at index k bounds cells k-1 and k

Cell-Centered Fields (nvert):

• Shifted indexing: Physical cell c stored at nvert[c+1]
• For mx=26 (25 cells): Cell 0→nvert[1], Cell 23→nvert[24]
• For mz=98 (96 cells): Cell 0→nvert[1], Cell 95→nvert[96]
• nvert[0] and nvert[mx-1] are ghost values

Face-to-Index Mapping

Example for a domain with mx=26, my=26, mz=98:

Face ID	Node Index	Face Metric	Adjacent Cell (shifted)	Physical Extent
BC_FACE_NEG_X	i=0	csi[k][j][0]	nvert[k][j][1] (Cell 0)	j∈[0,24], k∈[0,96]
BC_FACE_POS_X	i=24	csi[k][j][24]	nvert[k][j][24] (Cell 23)	j∈[0,24], k∈[0,96]
BC_FACE_NEG_Y	j=0	eta[k][0][i]	nvert[k][1][i] (Cell 0)	i∈[0,24], k∈[0,96]
BC_FACE_POS_Y	j=24	eta[k][24][i]	nvert[k][24][i] (Cell 23)	i∈[0,24], k∈[0,96]
BC_FACE_NEG_Z	k=0	zet[0][j][i]	nvert[1][j][i] (Cell 0)	i∈[0,24], j∈[0,24]
BC_FACE_POS_Z	k=96	zet[96][j][i]	nvert[96][j][i] (Cell 95)	i∈[0,24], j∈[0,24]

Algorithm

The function performs two separate computations with different loop bounds:

1. Center Calculation (uses ALL physical nodes):

• Loop over all physical nodes on the face (excluding dummy indices)
• Accumulate coordinate sums: Σx, Σy, Σz
• Count number of nodes
• Average: center = (Σx/n, Σy/n, Σz/n)

2. Area Calculation (uses INTERIOR cells only):

• Loop over interior cell range to avoid accessing ghost values in nvert
• For each face adjacent to a fluid cell (nvert < 0.1):
  • Compute area magnitude: |csi/eta/zet| = √(x² + y² + z²)
  • Accumulate to total area

Loop Bound Details

Why different bounds for center vs. area?

For BC_FACE_NEG_X at i=0 with my=26, mz=98:

Center calculation (coordinates):

• j ∈ [ys, j_max): Includes j=[0,24] (25 nodes), excludes dummy at j=25
• k ∈ [zs, k_max): Includes k=[0,96] (97 nodes), excludes dummy at k=97
• Total: 25 × 97 = 2,425 nodes

Area calculation (nvert checks):

• j ∈ [lys, lye): j=[1,24] (24 values), excludes boundaries
• k ∈ [lzs, lze): k=[1,96] (96 values), excludes boundaries
• Why restricted?
  • At j=0: nvert[k][0][1] is ghost (no cell at j=-1)
  • At j=25: nvert[k][25][1] is ghost (no cell at j=24, index 25 is dummy)
  • At k=0: nvert[0][j][1] is ghost (no cell at k=-1)
  • At k=97: nvert[97][j][1] is ghost (no cell at k=96, index 97 is dummy)
• Total: 24 × 96 = 2,304 interior cells adjacent to face

Area Calculation Formulas

Face area contributions are computed from metric tensor magnitudes:

• i-faces (±Xi): Area = |csi| = √(csi_x² + csi_y² + csi_z²)
• j-faces (±Eta): Area = |eta| = √(eta_x² + eta_y² + eta_z²)
• k-faces (±Zeta): Area = |zet| = √(zet_x² + zet_y² + zet_z²)

Parameters

[in]	user	Pointer to UserCtx containing grid info, DMs, and field vectors
[in]	face_id	Enum identifying which boundary face to analyze (BC_FACE_NEG_X, etc.)
[out]	face_center	Pointer to Cmpnts structure to store computed geometric center (x,y,z)
[out]	face_area	Pointer to PetscReal to store computed total face area

Returns

PetscErrorCode Returns 0 on success, non-zero PETSc error code on failure

Note

This function uses MPI_Allreduce, so it must be called collectively by all ranks

Only ranks that own the specified boundary face contribute to the calculation

Center calculation includes ALL physical nodes on the face

Area calculation ONLY includes faces adjacent to fluid cells (nvert < 0.1)

Dummy/unused indices (e.g., k=97, j=25 for standard test case) are excluded

Warning

Assumes grid and field arrays have been properly initialized

Incorrect face_id values will result in zero contribution from all ranks

See also

CanRankServiceFace() for determining rank ownership of boundary faces

BCFace enum for valid face_id values

LOG_FIELD_ANATOMY() for debugging field indexing

Example Usage:

Cmpnts inlet_center;
PetscReal inlet_area;
ierr = CalculateFaceCenterAndArea(user, BC_FACE_NEG_Z, &inlet_center, &inlet_area);
PetscPrintf(PETSC_COMM_WORLD, "Inlet center: (%.4f, %.4f, %.4f), Area: %.6f\n",
            inlet_center.x, inlet_center.y, inlet_center.z, inlet_area);

< Local sum of (x,y,z) coordinates

< Local sum of face area magnitudes

< Local count of nodes

< Global sum of coordinates

< Global sum of areas

< Global count of nodes

< i-range: [xs, xe)

< j-range: [ys, ye)

< k-range: [zs, ze)

< Physical domain size in i (exclude dummy)

< Physical domain size in j (exclude dummy)

< Physical domain size in k (exclude dummy)

< Start at 1 if on -Xi boundary

< End at mx-1 if on +Xi boundary

< Start at 1 if on -Eta boundary

< End at my-1 if on +Eta boundary

< Start at 1 if on -Zeta boundary

< End at mz-1 if on +Zeta boundary

< Exclude dummy at i=mx-1 (e.g., i=25)

< Exclude dummy at j=my-1 (e.g., j=25)

< Exclude dummy at k=mz-1 (e.g., k=97)

< Local ghosted coordinate vector

< Nodal coordinates [k][j][i]

< Face metric tensors [k][j][i]

< Cell blanking field [k][j][i] (shifted +1)

Definition at line 1276 of file grid.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    
    // ========================================================================
    // Local accumulators for this rank's contribution
    // ========================================================================
    PetscReal      local_sum[3] = {0.0, 0.0, 0.0};  ///< Local sum of (x,y,z) coordinates
    PetscReal      localAreaSum = 0.0;               ///< Local sum of face area magnitudes
    PetscCount     local_n_points = 0;               ///< Local count of nodes
    
    // ========================================================================
    // Global accumulators after MPI reduction
    // ========================================================================
    PetscReal      global_sum[3] = {0.0, 0.0, 0.0};  ///< Global sum of coordinates
    PetscReal      globalAreaSum = 0.0;              ///< Global sum of areas
    PetscCount     global_n_points = 0;              ///< Global count of nodes
    
    // ========================================================================
    // Grid information and array pointers
    // ========================================================================
    DM             da = user->da;
    info = user->info;
    
    // Rank's owned range in global indices
    PetscInt       xs = info.xs, xe = info.xs + info.xm;  ///< i-range: [xs, xe)
    PetscInt       ys = info.ys, ye = info.ys + info.ym;  ///< j-range: [ys, ye)
    PetscInt       zs = info.zs, ze = info.zs + info.zm;  ///< k-range: [zs, ze)
    
    // Global domain dimensions (total allocated, includes dummy at end)
    PetscInt       mx = info.mx, my = info.my, mz = info.mz;
    PetscInt       IM = user->IM;  ///< Physical domain size in i (exclude dummy)
    PetscInt       JM = user->JM;  ///< Physical domain size in j (exclude dummy)
    PetscInt       KM = user->KM;  ///< Physical domain size in k (exclude dummy)
 
    // ========================================================================
    // Interior loop bounds (adjusted to avoid ghost/boundary cells)
    // These are used for nvert checks where we need valid cell indices
    // ========================================================================
    PetscInt lxs = xs; if(xs == 0)  lxs = xs + 1;   ///< Start at 1 if on -Xi boundary
    PetscInt lxe = xe; if(xe == mx) lxe = xe - 1;   ///< End at mx-1 if on +Xi boundary
    PetscInt lys = ys; if(ys == 0)  lys = ys + 1;   ///< Start at 1 if on -Eta boundary
    PetscInt lye = ye; if(ye == my) lye = ye - 1;   ///< End at my-1 if on +Eta boundary
    PetscInt lzs = zs; if(zs == 0)  lzs = zs + 1;   ///< Start at 1 if on -Zeta boundary
    PetscInt lze = ze; if(ze == mz) lze = ze - 1;   ///< End at mz-1 if on +Zeta boundary
 
    // ========================================================================
    // Physical node bounds (exclude dummy indices at mx-1, my-1, mz-1)
    // These are used for coordinate loops where we want ALL physical nodes
    // ========================================================================
    PetscInt i_max = (xe == mx) ? mx - 1 : xe;  ///< Exclude dummy at i=mx-1 (e.g., i=25)
    PetscInt j_max = (ye == my) ? my - 1 : ye;  ///< Exclude dummy at j=my-1 (e.g., j=25)
    PetscInt k_max = (ze == mz) ? mz - 1 : ze;  ///< Exclude dummy at k=mz-1 (e.g., k=97)
 
    // ========================================================================
    // Array pointers for field access
    // ========================================================================
    Vec            lCoor;                       ///< Local ghosted coordinate vector
    Cmpnts         ***coor;                    ///< Nodal coordinates [k][j][i]
    Cmpnts         ***csi, ***eta, ***zet;    ///< Face metric tensors [k][j][i]
    PetscReal      ***nvert;                   ///< Cell blanking field [k][j][i] (shifted +1)
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // ========================================================================
    // Step 1: Check if this rank owns the specified boundary face
    // ========================================================================
    PetscBool owns_face = PETSC_FALSE;
    ierr = CanRankServiceFace(&info,IM,JM,KM,face_id,&owns_face); CHKERRQ(ierr);
    if(owns_face){    
        // ========================================================================
        // Step 2: Get read-only array access for all required fields
        // ========================================================================
        ierr = DMGetCoordinatesLocal(user->da, &lCoor); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, lCoor, &coor); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->da, user->lNvert, &nvert); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
 
        // ========================================================================
        // Step 3: Loop over the specified face and accumulate center and area
        // ========================================================================
        switch (face_id) {
            
            // ====================================================================
            // BC_FACE_NEG_X: Face at i=0 (bottom boundary in i-direction)
            // ====================================================================
            case BC_FACE_NEG_X:
                if (xs == 0) {
                    PetscInt i = 0;  // Face is at node index i=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // Loop over ALL physical nodes on this face
                    // For my=26, mz=98: j∈[0,24], k∈[0,96] → 25×97 = 2,425 nodes
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt j = ys; j < j_max; j++) {
                            // Accumulate coordinates at node [k][j][0]
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // Loop over interior range where nvert checks are valid
                    // For my=26, mz=98: j∈[1,24], k∈[1,96] → 24×96 = 2,304 cells
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt j = lys; j < lye; j++) {
                            // Check if adjacent cell is fluid
                            // nvert[k][j][i+1] = nvert[k][j][1] checks Cell 0
                            // (Physical Cell 0 in j-k plane, stored at shifted index [1])
                            if (nvert[k][j][i+1] < 0.1) {
                                // Cell is fluid - add face area contribution
                                // Face area = magnitude of csi metric at [k][j][0]
                                localAreaSum += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                                csi[k][j][i].y * csi[k][j][i].y +
                                                csi[k][j][i].z * csi[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_X: Face at i=IM-1 (top boundary in i-direction)
            // ====================================================================
            case BC_FACE_POS_X:
                if (xe == mx) {
                    PetscInt i = mx - 2;  // Last physical node (e.g., i=24 for mx=26)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt j = ys; j < j_max; j++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt j = lys; j < lye; j++) {
                            // Check if adjacent cell is fluid
                            // nvert[k][j][i] = nvert[k][j][24] checks last cell (Cell 23)
                            // (Physical Cell 23, stored at shifted index [24])
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of csi metric at [k][j][24]
                                localAreaSum += sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                                csi[k][j][i].y * csi[k][j][i].y +
                                                csi[k][j][i].z * csi[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_NEG_Y: Face at j=0 (bottom boundary in j-direction)
            // ====================================================================
            case BC_FACE_NEG_Y:
                if (ys == 0) {
                    PetscInt j = 0;  // Face is at node index j=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, mz=98: i∈[0,24], k∈[0,96] → 25×97 = 2,425 nodes
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, mz=98: i∈[1,24], k∈[1,96] → 24×96 = 2,304 cells
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j+1][i] = nvert[k][1][i] checks Cell 0
                            if (nvert[k][j+1][i] < 0.1) {
                                // Face area = magnitude of eta metric at [k][0][i]
                                localAreaSum += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                                eta[k][j][i].y * eta[k][j][i].y +
                                                eta[k][j][i].z * eta[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_Y: Face at j=JM-1 (top boundary in j-direction)
            // ====================================================================
            case BC_FACE_POS_Y:
                if (ye == my) {
                    PetscInt j = my - 2;  // Last physical node (e.g., j=24 for my=26)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    for (PetscInt k = zs; k < k_max; k++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    for (PetscInt k = lzs; k < lze; k++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j][i] = nvert[k][24][i] checks last cell (Cell 23)
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of eta metric at [k][24][i]
                                localAreaSum += sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                                eta[k][j][i].y * eta[k][j][i].y +
                                                eta[k][j][i].z * eta[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_NEG_Z: Face at k=0 (inlet, bottom boundary in k-direction)
            // ====================================================================
            case BC_FACE_NEG_Z:
                if (zs == 0) {
                    PetscInt k = 0;  // Face is at node index k=0
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, my=26: i∈[0,24], j∈[0,24] → 25×25 = 625 nodes
                    for (PetscInt j = ys; j < j_max; j++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, my=26: i∈[1,24], j∈[1,24] → 24×24 = 576 cells
                    for (PetscInt j = lys; j < lye; j++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k+1][j][i] = nvert[1][j][i] checks Cell 0
                            // (Physical Cell 0 in i-j plane, stored at shifted index [1])
                            if (nvert[k+1][j][i] < 0.1) {
                                // Face area = magnitude of zet metric at [0][j][i]
                                localAreaSum += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                                zet[k][j][i].y * zet[k][j][i].y +
                                                zet[k][j][i].z * zet[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            // ====================================================================
            // BC_FACE_POS_Z: Face at k=KM-1 (outlet, top boundary in k-direction)
            // ====================================================================
            case BC_FACE_POS_Z:
                if (ze == mz) {
                    PetscInt k = mz - 2;  // Last physical node (e.g., k=96 for mz=98)
                    
                    // ---- Part 1: Center calculation (ALL physical nodes) ----
                    // For mx=26, my=26: i∈[0,24], j∈[0,24] → 25×25 = 625 nodes
                    for (PetscInt j = ys; j < j_max; j++) {
                        for (PetscInt i = xs; i < i_max; i++) {
                            local_sum[0] += coor[k][j][i].x;
                            local_sum[1] += coor[k][j][i].y;
                            local_sum[2] += coor[k][j][i].z;
                            local_n_points++;
                        }
                    }
                    
                    // ---- Part 2: Area calculation (INTERIOR cells only) ----
                    // For mx=26, my=26: i∈[1,24], j∈[1,24] → 24×24 = 576 cells
                    for (PetscInt j = lys; j < lye; j++) {
                        for (PetscInt i = lxs; i < lxe; i++) {
                            // nvert[k][j][i] = nvert[96][j][i] checks last cell (Cell 95)
                            // (Physical Cell 95, stored at shifted index [96])
                            if (nvert[k][j][i] < 0.1) {
                                // Face area = magnitude of zet metric at [96][j][i]
                                localAreaSum += sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                                zet[k][j][i].y * zet[k][j][i].y +
                                                zet[k][j][i].z * zet[k][j][i].z);
                            }
                        }
                    }
                }
                break;
                
            default:
                SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 
                        "Unknown face_id %d in CalculateFaceCenterAndArea", face_id);
        }
        
        // ========================================================================
        // Step 4: Restore array access (release pointers)
        // ========================================================================
        ierr = DMDAVecRestoreArrayRead(user->fda, lCoor, &coor); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, &nvert); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    }    
    // ========================================================================
    // Step 5: Perform MPI reductions to get global sums
    // ========================================================================
    // Sum coordinate contributions from all ranks
    ierr = MPI_Allreduce(local_sum, global_sum, 3, MPI_DOUBLE, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // Sum node counts from all ranks
    ierr = MPI_Allreduce(&local_n_points, &global_n_points, 1, MPI_COUNT, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    
    // Sum area contributions from all ranks
    ierr = MPI_Allreduce(&localAreaSum, &globalAreaSum, 1, MPI_DOUBLE, MPI_SUM, 
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    // ========================================================================
    // Step 6: Calculate geometric center by averaging coordinates
    // ========================================================================
    if (global_n_points > 0) {
        face_center->x = global_sum[0] / global_n_points;
        face_center->y = global_sum[1] / global_n_points;
        face_center->z = global_sum[2] / global_n_points;
        LOG_ALLOW(GLOBAL, LOG_INFO, 
                  "Calculated center for Face %s: (x=%.4f, y=%.4f, z=%.4f) from %lld nodes\n", 
                  BCFaceToString(face_id), 
                  face_center->x, face_center->y, face_center->z,
                  (long long)global_n_points);
    } else {
        // No nodes found - this should not happen for a valid face
        LOG_ALLOW(GLOBAL, LOG_WARNING, 
                  "WARNING: Face %s identified but no grid points found. Center not calculated.\n",
                  BCFaceToString(face_id));
        face_center->x = face_center->y = face_center->z = 0.0;
    }
    
    // ========================================================================
    // Step 7: Return computed total area
    // ========================================================================
    *face_area = globalAreaSum;
    LOG_ALLOW(GLOBAL, LOG_INFO, 
              "Calculated area for Face %s: Area=%.6f\n", 
              BCFaceToString(face_id), *face_area);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: