Internal Scattering Helpers

Lower-level functions used by the main scattering routines. More...

Collaboration diagram for Internal Scattering Helpers:

Macros
#define	__FUNCT__   "InterpolateFieldFromCornerToCenter_Vector"
#define	__FUNCT__   "InterpolateFieldFromCornerToCenter_Scalar"
#define	__FUNCT__   "TestCornerToCenterInterpolation"
#define	__FUNCT__   "InterpolateFieldFromCenterToCorner_Vector"
#define	__FUNCT__   "InterpolateFieldFromCenterToCorner_Scalar"
#define	__FUNCT__   "InterpolateEulerFieldFromCenterToSwarm"
#define	__FUNCT__   "InterpolateEulerFieldFromCornerToSwarm"
#define	__FUNCT__   "InterpolateEulerFieldToSwarm"
#define	__FUNCT__   "InterpolateAllFieldsToSwarm"
#define	__FUNCT__   "AccumulateParticleField"
#define	__FUNCT__   "NormalizeGridVectorByCount"
#define	__FUNCT__   "ScatterParticleFieldToEulerField_Internal"
#define	__FUNCT__   "ScatterParticleFieldToEulerField"
#define	__FUNCT__   "ScatterAllParticleFieldsToEulerFields"
#define	__FUNCT__   "InterpolateCornerToFaceCenter_Scalar"
#define	__FUNCT__   "InterpolateCornerToFaceCenter_Vector"

Functions
PetscErrorCode	AccumulateParticleField (DM swarm, ParticleFieldId particle_field_id, DM gridSumDM, Vec gridSumVec)
	Accumulates a particle field (scalar or vector) into a target grid sum vector.
PetscErrorCode	NormalizeGridVectorByCount (DM countDM, Vec countVec, DM dataDM, Vec sumVec, Vec avgVec)
	Normalizes a grid vector of sums by a grid vector of counts to produce an average.

Detailed Description

Lower-level functions used by the main scattering routines.

Macro Definition Documentation

__FUNCT__ [1/16]

#define __FUNCT__   "InterpolateFieldFromCornerToCenter_Vector"

Definition at line 20 of file interpolation.c.

__FUNCT__ [2/16]

#define __FUNCT__   "InterpolateFieldFromCornerToCenter_Scalar"

Definition at line 20 of file interpolation.c.

__FUNCT__ [3/16]

#define __FUNCT__   "TestCornerToCenterInterpolation"

Definition at line 20 of file interpolation.c.

__FUNCT__ [4/16]

#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Vector"

Definition at line 20 of file interpolation.c.

__FUNCT__ [5/16]

#define __FUNCT__   "InterpolateFieldFromCenterToCorner_Scalar"

Definition at line 20 of file interpolation.c.

__FUNCT__ [6/16]

#define __FUNCT__   "InterpolateEulerFieldFromCenterToSwarm"

Definition at line 20 of file interpolation.c.

__FUNCT__ [7/16]

#define __FUNCT__   "InterpolateEulerFieldFromCornerToSwarm"

Definition at line 20 of file interpolation.c.

__FUNCT__ [8/16]

#define __FUNCT__   "InterpolateEulerFieldToSwarm"

Definition at line 20 of file interpolation.c.

__FUNCT__ [9/16]

#define __FUNCT__   "InterpolateAllFieldsToSwarm"

Definition at line 20 of file interpolation.c.

__FUNCT__ [10/16]

#define __FUNCT__   "AccumulateParticleField"

Definition at line 20 of file interpolation.c.

__FUNCT__ [11/16]

#define __FUNCT__   "NormalizeGridVectorByCount"

Definition at line 20 of file interpolation.c.

__FUNCT__ [12/16]

#define __FUNCT__   "ScatterParticleFieldToEulerField_Internal"

Definition at line 20 of file interpolation.c.

__FUNCT__ [13/16]

#define __FUNCT__   "ScatterParticleFieldToEulerField"

Definition at line 20 of file interpolation.c.

__FUNCT__ [14/16]

#define __FUNCT__   "ScatterAllParticleFieldsToEulerFields"

Definition at line 20 of file interpolation.c.

__FUNCT__ [15/16]

#define __FUNCT__   "InterpolateCornerToFaceCenter_Scalar"

Definition at line 20 of file interpolation.c.

__FUNCT__ [16/16]

#define __FUNCT__   "InterpolateCornerToFaceCenter_Vector"

Definition at line 20 of file interpolation.c.

Function Documentation

AccumulateParticleField()

PetscErrorCode AccumulateParticleField	(	DM	swarm,
		ParticleFieldId	particle_field_id,
		DM	gridSumDM,
		Vec	gridSumVec
	)

Accumulates a particle field (scalar or vector) into a target grid sum vector.

This function iterates through local particles, identifies their cell using the "DMSwarm_CellID" field, and adds the particle's field value (particleFieldName) to the corresponding cell location in the gridSumVec. It handles both scalar (DOF=1) and vector (DOF=3) fields automatically based on the DOF of gridSumDM.

IMPORTANT: The caller must ensure gridSumVec is zeroed before calling this function if a fresh sum calculation is desired.

Parameters

[in]	swarm	The DMSwarm containing particles.
[in]	particle_field_id	Typed identity of the particle field (must match DOF).
[in]	gridSumDM	The DMDA associated with gridSumVec. Its DOF determines how many components are accumulated.
[in,out]	gridSumVec	The Vec (associated with gridSumDM) to accumulate sums into.

Returns

PetscErrorCode 0 on success. Errors if fields don't exist or DMs are incompatible.

Definition at line 1478 of file interpolation.c.

{
    PetscErrorCode    ierr;
    PetscInt          dof;
    PetscInt          nlocal, p;
    const PetscReal   *particle_arr = NULL;
    const PetscInt    *cell_id_arr = NULL;
    
    // DMDA Accessors
    PetscScalar       ***arr_1d = NULL;      // For scalar fields
    PetscScalar       ****arr_3d = NULL;     // For vector fields
    
    // Ghosted dimension variables
    PetscInt          gxs, gys, gzs, gxm, gym, gzm;
    PetscMPIInt       rank;
    char              msg[ERROR_MSG_BUFFER_SIZE];
    const ParticleFieldDescriptor *descriptor = NULL;
    const char        *particleFieldName = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    ierr = ParticleFieldGetDescriptor(particle_field_id, &descriptor); CHKERRQ(ierr);
    PetscCheck(descriptor->data_type == PETSC_REAL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
               "Particle-to-grid accumulation requires PETSC_REAL data; field '%s' uses %s.",
               descriptor->canonical_name, PetscDataTypes[descriptor->data_type]);
    particleFieldName = descriptor->canonical_name;
 
    // --- 1. Validation & Setup ---
    if (!swarm || !gridSumDM || !localAccumulatorVec) 
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null input in AccumulateParticleField.");
 
    // Get DMDA information
    ierr = DMDAGetInfo(gridSumDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    // Get Ghosted Corners (Global indices of the ghosted patch start, and its dimensions)
    ierr = DMDAGetGhostCorners(gridSumDM, &gxs, &gys, &gzs, &gxm, &gym, &gzm); CHKERRQ(ierr);
 
    // --- 2. Verify Vector Type (Global vs Local) ---
    {
        PetscInt vecSize;
        PetscInt expectedLocalSize = gxm * gym * gzm * dof;
        ierr = VecGetSize(localAccumulatorVec, &vecSize); CHKERRQ(ierr);
        
        if (vecSize != expectedLocalSize) {
            PetscSNPrintf(msg, sizeof(msg), 
                "Vector dimension mismatch! Expected Ghosted Local Vector size %d (gxm*gym*gzm*dof), got %d. "
                "Did you pass a Global Vector instead of a Local Vector?", 
                expectedLocalSize, vecSize);
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg);
        }
    }
 
    // --- 3. Acquire Particle Data ---
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    // These calls will fail nicely if the field doesn't exist
    ierr = DMSwarmGetField(swarm, particleFieldName, NULL, NULL, (void **)&particle_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void **)&cell_id_arr); CHKERRQ(ierr);
 
    // --- 4. Acquire Grid Accessors ---
    // DMDAVecGetArray* handles the mapping from Global (i,j,k) to the underlying Local Array index
    if (dof == 1) {
        ierr = DMDAVecGetArray(gridSumDM, localAccumulatorVec, &arr_1d); CHKERRQ(ierr);
    } else if (dof == 3) {
        ierr = DMDAVecGetArrayDOF(gridSumDM, localAccumulatorVec, &arr_3d); CHKERRQ(ierr);
    } else {
         PetscSNPrintf(msg, sizeof(msg), "Unsupported DOF=%d. AccumulateParticleField supports DOF 1 or 3.", dof);
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "%s", msg);
    }
 
    // --- 5. Accumulate Loop ---
    // Iterate through all local particles
    for (p = 0; p < nlocal; ++p) {
        // Retrieve Geometric Grid Index (0-based)
        PetscInt i_geom = cell_id_arr[p * 3 + 0];
        PetscInt j_geom = cell_id_arr[p * 3 + 1];
        PetscInt k_geom = cell_id_arr[p * 3 + 2];
 
        // Apply Shift (+1) to match Memory Layout (Index 0 is boundary/ghost)
        PetscInt i = i_geom + 1;
        PetscInt j = j_geom + 1;
        PetscInt k = k_geom + 1;
 
        // Bounds Check: Ensure (i,j,k) falls within the Local Ghosted Patch.
        // This allows writing to ghost slots which will later be reduced to the owner rank.
        if (i >= gxs && i < gxs + gxm && 
            j >= gys && j < gys + gym && 
            k >= gzs && k < gzs + gzm)
        {
            if (dof == 1) {
                arr_1d[k][j][i] += particle_arr[p];
            } else {
                // For DOF=3, unroll the loop for slight optimization
                arr_3d[k][j][i][0] += particle_arr[p * 3 + 0];
                arr_3d[k][j][i][1] += particle_arr[p * 3 + 1];
                arr_3d[k][j][i][2] += particle_arr[p * 3 + 2];
            }
        } 
        // Note: Particles outside the ghost layer are skipped. This is expected behavior
        // if particles have not yet been migrated or localized correctly.
    }
 
    // --- 6. Restore Arrays and Fields ---
    if (dof == 1) {
        ierr = DMDAVecRestoreArray(gridSumDM, localAccumulatorVec, &arr_1d); CHKERRQ(ierr);
    } else {
        ierr = DMDAVecRestoreArrayDOF(gridSumDM, localAccumulatorVec, &arr_3d); CHKERRQ(ierr);
    }
    
    ierr = DMSwarmRestoreField(swarm, particleFieldName, NULL, NULL, (void **)&particle_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void **)&cell_id_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

NormalizeGridVectorByCount()

PetscErrorCode NormalizeGridVectorByCount	(	DM	countDM,
		Vec	countVec,
		DM	dataDM,
		Vec	sumVec,
		Vec	avgVec
	)

Normalizes a grid vector of sums by a grid vector of counts to produce an average.

Calculates avgVec[i] = sumVec[i] / countVec[i] for each component of each OWNED cell where countVec[i] > 0. Sets avgVec[i] = 0 otherwise. Handles both scalar (DOF=1) and vector (DOF=3) data fields based on dataDM. Uses basic VecGetArray/VecGetArrayRead and manual index calculation.

Parameters

[in]	countDM	The DMDA associated with countVec (must have DOF=1).
[in]	countVec	The Vec containing particle counts per cell (read-only).
[in]	dataDM	The DMDA associated with sumVec and avgVec (must have DOF=1 or DOF=3).
[in]	sumVec	The Vec containing the accumulated sums per cell (read-only).
[in,out]	avgVec	The Vec where the calculated averages will be stored (overwritten).

Returns

PetscErrorCode 0 on success.

Definition at line 1600 of file interpolation.c.

{
    PetscErrorCode ierr;
    PetscInt       data_dof;
    PetscInt       count_dof;
    PetscMPIInt    rank;
    char           msg[ERROR_MSG_BUFFER_SIZE];
 
    // Pointers for DMDA array accessors - declare specific types
    PetscScalar    ***count_arr_3d = NULL;      // For DOF=1 count vector (3D DMDA)
    PetscScalar    ***sum_arr_scalar = NULL;    // For DOF=1 sum vector (3D DMDA)
    PetscScalar    ***avg_arr_scalar = NULL;    // For DOF=1 avg vector (3D DMDA)
    PetscScalar    ****sum_arr_vector = NULL;   // For DOF=3 sum vector (3D DMDA + DOF)
    PetscScalar    ****avg_arr_vector = NULL;   // For DOF=3 avg vector (3D DMDA + DOF)
 
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- Validation ---
    ierr = DMDAGetInfo(countDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &count_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    ierr = DMDAGetInfo(dataDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &data_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, "%s", msg); }
    if (data_dof != 1 && data_dof != 3) { PetscSNPrintf(msg, sizeof(msg), "dataDM DOF must be 1 or 3, got %d.", data_dof); SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg); }
 
    // --- Get Array Access using appropriate DMDA accessors ---
    ierr = DMDAVecGetArrayRead(countDM, countVec, &count_arr_3d); CHKERRQ(ierr);
 
    if (data_dof == 1) {
        ierr = DMDAVecGetArrayRead(dataDM, sumVec, &sum_arr_scalar); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dataDM, avgVec, &avg_arr_scalar); CHKERRQ(ierr);
    } else { // data_dof == 3
        ierr = DMDAVecGetArrayDOFRead(dataDM, sumVec, &sum_arr_vector); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayDOF(dataDM, avgVec, &avg_arr_vector); CHKERRQ(ierr);
    }
 
    // Get the corners (global start indices) and dimensions of the *local owned* region
    PetscInt xs, ys, zs, xm, ym, zm;
    ierr = DMDAGetCorners(countDM, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ(ierr);
 
    // --- Normalize Over Owned Cells ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "(Rank %d): Normalizing DOF=%d data over owned range [%d:%d, %d:%d, %d:%d].\n",
              rank, data_dof, xs, xs+xm, ys, ys+ym, zs, zs+zm);
 
    // Loop using GLOBAL indices (i, j, k) over the range owned by this process
    for (PetscInt k = zs; k < zs + zm; ++k) {
        for (PetscInt j = ys; j < ys + ym; ++j) {
            for (PetscInt i = xs; i < xs + xm; ++i) {
 
                // Access the count using standard 3D indexing
                PetscScalar count = count_arr_3d[k][j][i];
 
                if (PetscRealPart(count) > 0.5) { // Use tolerance for float comparison
                    if (data_dof == 1) {
                        // Access scalar sum/avg using standard 3D indexing
                        avg_arr_scalar[k][j][i] = sum_arr_scalar[k][j][i] / count;
                    } else { // data_dof == 3
                        // Access vector components using DOF indexing on the last dimension
                        for (PetscInt c = 0; c < data_dof; ++c) {
                            avg_arr_vector[k][j][i][c] = sum_arr_vector[k][j][i][c] / count;
                        }
                    }
                } else { // count is zero or negative
                    // Set average to zero
                    if (data_dof == 1) {
                        avg_arr_scalar[k][j][i] = 0.0;
                    } else { // data_dof == 3
                        for (PetscInt c = 0; c < data_dof; ++c) {
                            avg_arr_vector[k][j][i][c] = 0.0;
                        }
                    }
                } // end if count > 0.5
            } // end i loop
        } // end j loop
    } // end k loop
 
    // --- Restore Arrays using appropriate functions ---
    ierr = DMDAVecRestoreArrayRead(countDM, countVec, &count_arr_3d); CHKERRQ(ierr);
    if (data_dof == 1) {
        ierr = DMDAVecRestoreArrayRead(dataDM, sumVec, &sum_arr_scalar); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dataDM, avgVec, &avg_arr_scalar); CHKERRQ(ierr);
    } else { // data_dof == 3
        ierr = DMDAVecRestoreArrayDOFRead(dataDM, sumVec, &sum_arr_vector); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayDOF(dataDM, avgVec, &avg_arr_vector); CHKERRQ(ierr);
    }
 
    // --- Assemble Final Average Vector ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling final average vector (DOF=%d).\n", data_dof);
    ierr = VecAssemblyBegin(avgVec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(avgVec); CHKERRQ(ierr);
 
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function: