grid.h File Reference

Public interface for grid, solver, and metric setup routines. More...

#include "variables.h"
#include "logging.h"
#include "io.h"
#include "setup.h"
#include "AnalyticalSolutions.h"

Include dependency graph for grid.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	DefineAllGridDimensions (SimCtx *simCtx)
	Orchestrates the parsing and setting of grid dimensions for all blocks.
PetscErrorCode	InitializeAllGridDMs (SimCtx *simCtx)
	Orchestrates the creation of DMDA objects for every block and multigrid level.
PetscErrorCode	AssignAllGridCoordinates (SimCtx *simCtx)
	Orchestrates the assignment of physical coordinates to all DMDA objects.
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 information collected on rank 0 to all other 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.

Detailed Description

Public interface for grid, solver, and metric setup routines.

Definition in file grid.h.

Function Documentation

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:

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:

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:

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 of the local grid points owned by the current MPI process and stores the computed bounding box in the provided structure.

Parameters

[in]	user	Pointer to the user-defined context containing grid information.
[out]	localBBox	Pointer to the BoundingBox structure to store the computed bounding box.

Returns

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

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.

This function computes the local bounding box on each process, then collects all local bounding boxes on the root process (rank 0) using MPI. The result is stored in an array of BoundingBox structures on rank 0.

Parameters

[in]	user	Pointer to the user-defined context containing grid information.
[out]	allBBoxes	Pointer to a pointer where the array of gathered bounding boxes will be stored on rank 0. The caller on rank 0 must free this array.

Returns

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

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 information collected on rank 0 to all other ranks.

This function assumes that GatherAllBoundingBoxes() was previously called, so bboxlist is allocated and populated on rank 0. All other ranks will allocate memory for bboxlist, and this function will use MPI_Bcast to distribute the bounding box data to them.

Parameters

[in]	user	Pointer to the UserCtx structure. (Currently unused in this function, but kept for consistency.)
[in,out]	bboxlist	Pointer to the array of BoundingBoxes. On rank 0, this should point to a valid array of size 'size' (where size is the number of MPI ranks). On non-root ranks, this function will allocate memory for bboxlist.

Returns

PetscErrorCode Returns 0 on success, non-zero on MPI or PETSc-related errors.

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: