Particle-to-Grid Scattering

Functions for scattering particle data (scalar or vector) onto Eulerian grid fields by averaging contributions within each cell. More...

Collaboration diagram for Particle-to-Grid Scattering:

Modules
	Internal Scattering Helpers
	Lower-level functions used by the main scattering routines.

Functions
PetscErrorCode	GetScatterTargetInfo (UserCtx *user, const char *particleFieldName, DM *targetDM, PetscInt *expected_dof)
	Determines the target Eulerian DM and expected DOF for scattering a given particle field.
PetscErrorCode	ScatterParticleFieldToEulerField (UserCtx *user, const char *particleFieldName, Vec eulerFieldAverageVec)
	Scatters a particle field (scalar or vector) to the corresponding Eulerian field average.
PetscErrorCode	ScatterAllParticleFieldsToEulerFields (UserCtx *user)
	Scatters a predefined set of particle fields to their corresponding Eulerian fields.
PetscErrorCode	CalculateParticleCountPerCell (UserCtx *user)
	Counts particles in each cell of the DMDA 'da' and stores the result in user->ParticleCount.

Detailed Description

Functions for scattering particle data (scalar or vector) onto Eulerian grid fields by averaging contributions within each cell.

This file provides a modular set of functions to perform particle-to-grid projection, specifically calculating cell-averaged quantities from particle properties. It assumes a PETSc environment using DMDA for the grids and DMSwarm for particles.

Key Features:

• Handles both scalar (DOF=1) and vector (DOF=3) particle fields.
• Uses a pre-calculated particle count vector (ParticleCount) for normalization.
• Implicitly determines the target Eulerian DM (typically user->da for scalars, user->fda for vectors) based on standard field names ("P", "Nvert", "Ucat", "Ucont","Psi").
• Modifies existing Eulerian field vectors (e.g., user->P, user->Ucat,user->Psi) in place.
• Provides a high-level wrapper function (ScatterAllParticleFieldsToEulerFields) to easily scatter a standard set of fields.
• Uses only the baseSETERRQ` macro for error reporting to maximize compiler compatibility.

Dependencies:

• PETSc library (DMDA, DMSwarm, Vec, IS, PetscLog, etc.)
• A UserCtx struct (defined elsewhere, e.g., "userctx.h") containing pointers to relevant DMs (da, fda), Vecs (ParticleCount, P, Nvert, Ucat, etc.), and the DMSwarm object (swarm).
• A custom DMSwarm field named "DMSwarm_CellID" (blockSize=3, type=PETSC_INT) must be registered and populated with the local cell indices for each particle.
• Logging infrastructure (LOG_ALLOW, etc.) assumed to be defined elsewhere.

Function Documentation

GetScatterTargetInfo()

PetscErrorCode GetScatterTargetInfo	(	UserCtx *	user,
		const char *	particleFieldName,
		DM *	targetDM,
		PetscInt *	expected_dof
	)

Determines the target Eulerian DM and expected DOF for scattering a given particle field.

Based on hardcoded rules mapping particle field names to user context DMs (da/fda). This function encapsulates the policy of where different fields should be scattered.

Parameters

[in]	user	Pointer to the UserCtx containing da and fda.
[in]	particleFieldName	Name of the particle field (e.g., "P", "Ucat").
[out]	targetDM	Pointer to store the determined target DM (da or fda).
[out]	expected_dof	Pointer to store the expected DOF (1 or 3) for this field.

Returns

PetscErrorCode Returns 0 on success, PETSC_ERR_ARG_UNKNOWN if field name is not recognized, or PETSC_ERR_ARG_NULL for NULL inputs.

Based on hardcoded rules mapping particle field names ("P", "Nvert", "Ucat", "Ucont") to user context DMs (user->da or user->fda). This function encapsulates the policy of where different fields should be scattered. Modify this function to add rules for custom fields.

Parameters

[in]	user	Pointer to the UserCtx containing required DMs (da, fda).
[in]	particleFieldName	Name of the particle field.
[out]	targetDM	Pointer to store the determined target DM (user->da or user->fda).
[out]	expected_dof	Pointer to store the expected DOF (1 or 3) for this field.

Returns

PetscErrorCode Returns 0 on success. Error codes:

• PETSC_ERR_ARG_NULL if required inputs are NULL.
• PETSC_ERR_ARG_WRONG if particleFieldName is not recognized.

Definition at line 1395 of file interpolation.c.

{
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- Input Validation ---
    // Check for NULL pointers in essential inputs
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL.");
    if (!user->da) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->da is NULL.");
    if (!user->fda) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->fda is NULL."); // Needed for vector fields
    if (!particleFieldName) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "particleFieldName is NULL.");
    if (!targetDM) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output targetDM pointer is NULL.");
    if (!expected_dof) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output expected_dof pointer is NULL.");
 
    // --- Determine Target DM and DOF based on Field Name ---
    // Compare the input field name with known scalar fields targeting 'da'
    if (strcmp(particleFieldName, "Psi") == 0 || strcmp(particleFieldName, "Nvert") == 0) {
        *expected_dof = 1;      // Scalar fields have DOF 1
        *targetDM = user->da;   // Target the primary scalar DMDA
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Field '%s' targets DM 'da' (DOF=1).\n", particleFieldName);
    }
    // Compare with known vector fields targeting 'fda'
    else if (strcmp(particleFieldName, "Ucat") == 0 || strcmp(particleFieldName, "Ucont") == 0) {
        *expected_dof = 3;      // Vector fields have DOF 3
        *targetDM = user->fda;  // Target the vector DMDA (often node-based)
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Field '%s' targets DM 'fda' (DOF=3).\n", particleFieldName);
    }
    // --- Add rules for other fields here ---
    // else if (strcmp(particleFieldName, "SomeOtherScalar") == 0) { *expected_dof = 1; *targetDM = user->da; }
    // else if (strcmp(particleFieldName, "SomeOtherVector") == 0) { *expected_dof = 3; *targetDM = user->someOtherDM; }
    else {
        // The provided field name doesn't match any known rules
        *targetDM = NULL; // Indicate failure
        *expected_dof = 0;
        // Format the error message manually
        PetscSNPrintf(msg, sizeof(msg), "Field name '%s' is not recognized for automatic DM selection.", particleFieldName);
        // Use SETERRQ with the formatted message and appropriate error code
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, msg); // Use WRONG argument error code
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ScatterParticleFieldToEulerField()

PetscErrorCode ScatterParticleFieldToEulerField	(	UserCtx *	user,
		const char *	particleFieldName,
		Vec	eulerFieldAverageVec
	)

Scatters a particle field (scalar or vector) to the corresponding Eulerian field average.

This is the main user-facing function. It determines the target Eulerian DM based on the particleFieldName, validates the provided eulerFieldAverageVec against the target DM, and then orchestrates the scatter operation by calling the internal helper function ScatterParticleFieldToEulerField_Internal. The final averaged result is stored IN-PLACE in eulerFieldAverageVec.

Parameters

[in]	user	Pointer to UserCtx containing da, fda, swarm, ParticleCount.
[in]	particleFieldName	Name of the field in the DMSwarm (e.g., "P", "Ucat").
[in,out]	eulerFieldAverageVec	Pre-created Vec associated with the correct target DM (implicitly da or fda). Result stored here.

Returns

PetscErrorCode 0 on success. Errors on NULL input, unrecognized field name, or incompatible target vector.

This is the primary user-facing function for scattering a single field. It determines the target Eulerian DM (user->da or user->fda) based on the particleFieldName, validates that the provided eulerFieldAverageVec is compatible with that target DM, and then orchestrates the scatter operation by calling an internal helper function. The final averaged result is stored IN-PLACE in the eulerFieldAverageVec.

Parameters

[in]	user	Pointer to UserCtx containing da, fda, swarm, ParticleCount.
[in]	particleFieldName	Name of the field in the DMSwarm (e.g., "Psi", "Ucat").
[in,out]	eulerFieldAverageVec	Pre-created Vec associated with the correct target DM (implicitly da or fda). Result stored here. This vector should be zeroed by the caller if desired before this function (e.g., using VecSet). The wrapper ScatterAllParticleFieldsToEulerFields handles this zeroing.

Returns

PetscErrorCode 0 on success. Errors on NULL input, unrecognized field name, incompatible target vector, or if the particle field doesn't exist.

Definition at line 1798 of file interpolation.c.

{
    PetscErrorCode ierr;
    DM             targetDM = NULL;       // Will point to user->da or user->fda
    PetscInt       expected_dof = 0;      // Will be 1 or 3
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- Essential Input Validation ---
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL.");
    if (!user->swarm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->swarm is NULL.");
    if (!user->ParticleCount) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->ParticleCount is NULL.");
    if (!eulerFieldAverageVec) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output eulerFieldAverageVec is NULL.");
    // particleFieldName validity checked within GetScatterTargetInfo
 
    // --- Determine Target DM & DOF using the helper function ---
    ierr = GetScatterTargetInfo(user, particleFieldName, &targetDM, &expected_dof); CHKERRQ(ierr);
    // If GetScatterTargetInfo failed (e.g., unknown name), CHKERRQ would have returned.
 
    // --- Validate the provided Target Vec's Compatibility ---
    DM vec_dm;
    PetscInt vec_dof;
    // Check that the provided average vector has a DM associated with it
    ierr = VecGetDM(eulerFieldAverageVec, &vec_dm); CHKERRQ(ierr);
    if (!vec_dm) {
         PetscSNPrintf(msg, sizeof(msg), "Provided eulerFieldAverageVec for field '%s' does not have an associated DM.", particleFieldName);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, msg);
    }
    // Get the block size (DOF) of the provided vector
    ierr = VecGetBlockSize(eulerFieldAverageVec, &vec_dof); CHKERRQ(ierr);
    // Compare the vector's associated DM with the one determined by the field name
    if (vec_dm != targetDM) {
        const char *target_dm_name = "targetDM", *vec_dm_name = "vec_dm";
        // Get actual names if possible for a more informative error message
        PetscObjectGetName((PetscObject)targetDM, &target_dm_name);
        PetscObjectGetName((PetscObject)vec_dm, &vec_dm_name);
        PetscSNPrintf(msg, sizeof(msg), "Provided eulerFieldAverageVec associated with DM '%s', but field '%s' requires scatter to DM '%s'.", vec_dm_name, particleFieldName, target_dm_name);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, msg);
    }
    // Compare the vector's DOF with the one expected for the field name
    if (vec_dof != expected_dof) {
         PetscSNPrintf(msg, sizeof(msg), "Field '%s' requires DOF %d, but provided eulerFieldAverageVec has DOF %d.", particleFieldName, expected_dof, vec_dof);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, msg);
    }
 
    // --- Perform Scatter using Internal Helper ---
    // Log intent before calling the core logic
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering field '%s' (DOF=%d).\n", particleFieldName, expected_dof);
    ierr = ScatterParticleFieldToEulerField_Internal(user, // Pass user context
                                                     particleFieldName, // Name of particle field
                                                     targetDM, // Determined target DM (da or fda)
                                                     expected_dof, // Determined DOF (1 or 3)
                                                     eulerFieldAverageVec); // The output vector
    CHKERRQ(ierr); // Handle potential errors from the internal function
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Successfully scattered field '%s'.\n", particleFieldName);
   
    PROFILE_FUNCTION_END;
   
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ScatterAllParticleFieldsToEulerFields()

PetscErrorCode ScatterAllParticleFieldsToEulerFields

(

UserCtx *

user

)

Scatters a predefined set of particle fields to their corresponding Eulerian fields.

This convenience function calls the unified ScatterParticleFieldToEulerField for a standard set of fields ("P", potentially others). It assumes the target Eulerian Vec objects (e.g., user->P, user->Ucat) exist in the UserCtx structure and are correctly associated with their respective DMs (user->da or user->fda). It zeros the target Vecs before scattering.

Parameters

[in,out]

user

Pointer to the UserCtx structure containing all required DMs, Vecs (ParticleCount, target Eulerian fields like P, Ucat), and swarm.

Returns

PetscErrorCode 0 on success. Errors if prerequisites (like ParticleCount) are missing or if underlying scatter calls fail.

This convenience function calls the unified ScatterParticleFieldToEulerField for a standard set of fields (currently just "Psi", others commented out). It assumes the target Eulerian Vec objects (e.g., user->Psi, user->Ucat) exist in the UserCtx structure and are correctly associated with their respective DMs (user->da or user->fda).

IMPORTANT: It zeros the target Vecs before scattering. Ensure user->ParticleCount has been computed accurately before calling this function, reflecting the current particle distribution.

Parameters

[in,out]

user

Pointer to the UserCtx structure containing all required DMs, Vecs (ParticleCount, target Eulerian fields like P, Ucat), and swarm.

Returns

PetscErrorCode 0 on success. Errors if prerequisites (like ParticleCount) are missing or if underlying scatter calls fail (e.g., particle field missing).

Definition at line 1887 of file interpolation.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting scattering of specified particle fields to Eulerian grids.\n");
 
    // --- Pre-computation Check: Ensure Particle Counts are Ready ---
    if (!user->ParticleCount) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "UserCtx->ParticleCount is NULL. Compute counts before calling ScatterAllParticleFieldsToEulerFields.");
    }
 
    // --- Scatter Particle Field "Psi" -> Eulerian Field user->Psi (on da) ---
    // Check if the target Eulerian vector 'user->Psi' exists.
    if (user->Psi) {
        
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Psi' to user->Psi.\n");
        // Zero the target vector before accumulating the new average for this step/call.
 
    // Debug Verification ------------------------------------------------
    Vec swarm_Psi;
    PetscReal Avg_Psi,Avg_swarm_Psi;
         
    ierr = VecMean(user->Psi,&Avg_Psi);
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Scalar(Psi) before scatter: %.4f.\n",Avg_Psi);
         
    ierr = DMSwarmCreateGlobalVectorFromField(user->swarm,"Psi",&swarm_Psi);
    ierr = VecMean(swarm_Psi,&Avg_swarm_Psi);
 
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Particle Scalar(Psi): %.4f.\n",Avg_swarm_Psi);
 
    ierr = DMSwarmDestroyGlobalVectorFromField(user->swarm,"Psi",&swarm_Psi);
    // Debug----------------------------------------------------------------
      
    //ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
        // Call the unified scatter function. It will handle DM determination and validation.
        // It will also error out if the *particle* field "Psi" doesn't exist in the swarm.
        ierr = ScatterParticleFieldToEulerField(user, "Psi", user->Psi); CHKERRQ(ierr);
    ierr = VecMean(user->Psi,&Avg_Psi);
    
    LOG_ALLOW(GLOBAL,LOG_DEBUG," Average of Scalar(Psi) after  scatter: %.4f.\n",Avg_Psi);
    } else {
        // Only log a warning if the target Eulerian field is missing in the context.
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Psi': UserCtx->Psi is NULL.\n");
    }
 
    // --- (Commented Out) Scatter Other Fields ---
    // To enable scattering for other fields, uncomment the relevant block
    // AND ensure:
    // 1. The corresponding particle field (e.g., "Nvert", "Ucat") exists in user->swarm.
    // 2. The corresponding target Eulerian Vec (e.g., user->Nvert, user->Ucat) exists in user->ctx.
    // 3. The target Eulerian Vec is associated with the correct DM (da for Nvert, fda for Ucat).
 
    /*
    // Scatter Particle Field "Nvert" -> Eulerian Field user->Nvert (on da)
    if (user->Nvert) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Nvert' to user->Nvert.\n");
        ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Nvert", user->Nvert); CHKERRQ(ierr);
    } else {
         LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Nvert': UserCtx->Nvert is NULL.\n");
    }
    */
 
    /*
    // Scatter Particle Field "Ucat" -> Eulerian Field user->Ucat (on fda)
    if (user->Ucat) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Ucat' to user->Ucat.\n");
        ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Ucat", user->Ucat); CHKERRQ(ierr);
    } else {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Ucat': UserCtx->Ucat is NULL.\n");
    }
    */
 
    /*
     // Scatter Particle Field "Ucont" -> Eulerian Field user->Ucont (on fda)
    if (user->Ucont) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scattering particle field 'Ucont' to user->Ucont.\n");
        ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
        ierr = ScatterParticleFieldToEulerField(user, "Ucont", user->Ucont); CHKERRQ(ierr);
    } else {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Skipping scatter for 'Ucont': UserCtx->Ucont is NULL.\n");
    }
    */
 
    // Add more fields as needed following the pattern above...
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished scattering specified particle fields.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

CalculateParticleCountPerCell()

PetscErrorCode CalculateParticleCountPerCell

(

UserCtx *

user

)

Counts particles in each cell of the DMDA 'da' and stores the result in user->ParticleCount.

Zeros the user->ParticleCount vector, then iterates through local particles. Reads the GLOBAL cell index (I, J, K) stored in the "DMSwarm_CellID" field. Uses DMDAVecGetArray to access the local portion of the count vector and increments the count at the global index (I, J, K) if it belongs to the local patch (including ghosts).

Parameters

[in,out]

user

Pointer to the UserCtx structure containing da, swarm, and ParticleCount.

Returns

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

Definition at line 485 of file ParticleMotion.c.

                                                            {
    PetscErrorCode ierr;
    DM             da = user->da;
    DM             swarm = user->swarm;
    Vec            countVec = user->ParticleCount;
    PetscInt       nlocal, p;
    PetscInt       *global_cell_id_arr; // Read GLOBAL cell IDs
    PetscScalar    ***count_arr_3d;     // Use 3D accessor
    PetscInt64       *PID_arr;
    PetscMPIInt    rank;
    char           msg[ERROR_MSG_BUFFER_SIZE];
    PetscInt       particles_counted_locally = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- Input Validation ---
    if (!da) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->da is NULL.");
    if (!swarm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->swarm is NULL.");
    if (!countVec) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx->ParticleCount is NULL.");
    // Check DOF of da
    PetscInt count_dof;
    ierr = DMDAGetInfo(da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &count_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
     if (count_dof != 1) { PetscSNPrintf(msg, sizeof(msg), "countDM must have DOF=1, got %d.", count_dof); SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, msg); }
 
    // --- Zero the count vector ---
    ierr = VecSet(countVec, 0.0); CHKERRQ(ierr);
 
    // --- Get Particle Data ---
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Accessing particle data.\n");
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm,"DMSwarm_CellID", NULL, NULL, (void **)&global_cell_id_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm,"DMSwarm_pid",NULL,NULL,(void **)&PID_arr);CHKERRQ(ierr);
 
    // --- Get Grid Vector Array using DMDA accessor ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Accessing ParticleCount vector array (using DMDAVecGetArray).\n");
    ierr = DMDAVecGetArray(da, countVec, &count_arr_3d); CHKERRQ(ierr);
 
    // Get local owned range for validation/logging if needed, but not for indexing with DMDAVecGetArray
    PetscInt xs, ys, zs, xm, ym, zm;
    ierr = DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ(ierr);
 
    // --- Accumulate Counts Locally ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "CalculateParticleCountPerCell (Rank %d): Processing %d local particles using GLOBAL CellIDs.\n",rank,nlocal);
    for (p = 0; p < nlocal; p++) {
        // Read the GLOBAL indices stored for this particle
        PetscInt64 i_geom = global_cell_id_arr[p * 3 + 0]; // Global i index
        PetscInt64 j_geom = global_cell_id_arr[p * 3 + 1]; // Global j index
        PetscInt64 k_geom = global_cell_id_arr[p * 3 + 2]; // Global k index
 
        // Apply the shift to ensure ParticleCount follows the indexing convention for cell-centered data in this codebase.
        PetscInt i = (PetscInt)i_geom + 1; // Shift for cell-centered
        PetscInt j = (PetscInt)j_geom + 1; // Shift for cell-centered
        PetscInt k = (PetscInt)k_geom + 1; // Shift for cell-centered
 
        // *** Bounds check is implicitly handled by DMDAVecGetArray for owned+ghost region ***
        // However, accessing outside this region using global indices WILL cause an error.
        // A preliminary check might still be wise if global IDs could be wild.
        // We rely on LocateAllParticles to provide valid global indices [0..IM-1] etc.
 
        LOG_LOOP_ALLOW(LOCAL,LOG_VERBOSE,p,100,"[Rank %d] Read CellID for p=%d, PID = %ld: (%ld, %ld, %ld)\n", rank, p,PID_arr[p],i, j, k);
    
        // *** Access the local array using GLOBAL indices ***
        // DMDAVecGetArray allows this, mapping global (I,J,K) to the correct
        // location within the local ghosted array segment.
        // This assumes (I,J,K) corresponds to a cell within the local owned+ghost region.
        // If (I,J,K) is for a cell owned by *another* rank (and not in the ghost region
        // of this rank), accessing count_arr_3d[K][J][I] will likely lead to a crash
        // or incorrect results. This highlights why storing LOCAL indices is preferred
        // for parallel runs. But proceeding with GLOBAL as requested:
        if (i >= xs && i < xs + xm  && 
            j >= ys && j < ys + ym  && // Adjust based on actual ghost width
            k >= zs && k < zs + zm  )   // This check prevents definite crashes but doesn't guarantee ownership
        {
 
             // Increment count at the location corresponding to GLOBAL index (I,J,K)
      //  LOG_ALLOW(LOCAL, LOG_DEBUG, "CalculateParticleCountPerCell (Rank %d): Particle %d with global CellID (%d, %d, %d) incremented with a particle.\n",rank, p, i, j, k);
             count_arr_3d[k][j][i] += 1.0;
             particles_counted_locally++;
         } else {
              // This particle's global ID is likely outside the range this rank handles (even ghosts)
              LOG_ALLOW(LOCAL, LOG_DEBUG, "CalculateParticleCountPerCell (Rank %d): Skipping particle %ld with global CellID (%ld, %ld, %ld) - likely outside local+ghost range.\n",rank, PID_arr[p] , i, j, k);
    }
    }
    LOG_ALLOW(LOCAL, LOG_DEBUG, "CalculateParticleCountPerCell (Rank %d): Local counting finished. Processed %d particles locally.\n", rank, particles_counted_locally);
 
    // --- Restore Access ---
    ierr = DMDAVecRestoreArray(da, countVec, &count_arr_3d); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm,"DMSwarm_CellID", NULL, NULL, (void **)&global_cell_id_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm,"DMSwarm_pid",NULL,NULL,(void **)&PID_arr);CHKERRQ(ierr);
 
    // --- Assemble Global Vector ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global ParticleCount vector.\n");
    ierr = VecAssemblyBegin(countVec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(countVec); CHKERRQ(ierr);
 
    // --- Verification Logging ---
    PetscReal total_counted_particles = 0.0, max_count_in_cell = 0.0;
    ierr = VecSum(countVec, &total_counted_particles); CHKERRQ(ierr);
    PetscInt max_idx_global = -1;
    ierr = VecMax(countVec, &max_idx_global, &max_count_in_cell); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Total counted globally = %.0f, Max count in cell = %.0f\n",
              total_counted_particles, max_count_in_cell);
 
    // --- ADD THIS DEBUGGING BLOCK ---
    if (max_idx_global >= 0) { // Check if VecMax found a location
         // Need to convert the flat global index back to 3D global index (I, J, K)
         // Get global grid dimensions (Nodes, NOT Cells IM/JM/KM)
         PetscInt M, N, P;
         ierr = DMDAGetInfo(da, NULL, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
         // Note: Assuming DOF=1 for countVec, index mapping uses node dimensions M,N,P from DMDA creation (IM+1, etc)
         // Re-check if your DMDA uses cell counts (IM) or node counts (IM+1) for Vec layout. Let's assume Node counts M,N,P.
         PetscInt Kmax = max_idx_global / (M * N);
         PetscInt Jmax = (max_idx_global % (M * N)) / M;
         PetscInt Imax = max_idx_global % M;
         LOG_ALLOW(GLOBAL, LOG_INFO, "  -> Max count located at global index (I,J,K) = (%d, %d, %d) [Flat index: %d]\n",
                   (int)Imax, (int)Jmax, (int)Kmax, (int)max_idx_global);
 
        // Also, let's explicitly check the count at (0,0,0)
        PetscScalar count_at_origin = 0.0;
        PetscScalar ***count_arr_for_check;
        ierr = DMDAVecGetArrayRead(da, countVec, &count_arr_for_check); CHKERRQ(ierr);
        // Check bounds before accessing - crucial if using global indices
        PetscInt xs, ys, zs, xm, ym, zm;
        ierr = DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ(ierr);
        if (0 >= xs && 0 < xs+xm && 0 >= ys && 0 < ys+ym && 0 >= zs && 0 < zs+zm) {
             count_at_origin = count_arr_for_check[0][0][0]; // Access using global index (0,0,0)
        } else {
            // Origin is not on this rank (relevant for parallel, but check anyway)
            count_at_origin = -999.0; // Indicate it wasn't accessible locally
        }
        ierr = DMDAVecRestoreArrayRead(da, countVec, &count_arr_for_check); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO, "  -> Count at global index (0,0,0) = %.1f\n", count_at_origin);
 
    } else {
         LOG_ALLOW(GLOBAL, LOG_WARNING, "  -> VecMax did not return a location for the maximum value.\n");
    }
    // --- END DEBUGGING BLOCK ---
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle counting complete.\n");
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function: