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	ScatterParticleFieldToEulerField (UserCtx *user, ParticleFieldId particle_field_id, 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.
static PetscErrorCode	ScatterParticleFieldToEulerField_Internal (UserCtx *user, ParticleFieldId particle_field_id, DM targetDM, PetscInt expected_dof, Vec eulerFieldAverageVec)
	Accumulate one particle field onto the Eulerian grid using the selected scatter stencil.

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.
• Uses particle-catalog metadata to resolve each supported Eulerian target.
• Modifies an existing, explicitly provided Eulerian field vector in place.
• Provides a high-level wrapper function (ScatterAllParticleFieldsToEulerFields) to easily scatter a standard set of fields.
• Uses only the base SETERRQ 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

ScatterParticleFieldToEulerField()

PetscErrorCode ScatterParticleFieldToEulerField	(	UserCtx *	user,
		ParticleFieldId	particle_field_id,
		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 from the particle catalog, 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]	particle_field_id	Typed identity of the DMSwarm field.
[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.

Definition at line 1802 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
    const ParticleFieldDescriptor *particle_descriptor = NULL;
    const FieldDescriptor *eulerian_descriptor = NULL;
    FieldView       target_view;
    const char     *particleFieldName = NULL;
 
    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.");
    ierr = ParticleFieldGetDescriptor(particle_field_id, &particle_descriptor); CHKERRQ(ierr);
    PetscCheck((particle_descriptor->capabilities & PARTICLE_FIELD_CAPABILITY_EULERIAN_SCATTER) != 0,
               PETSC_COMM_SELF, PETSC_ERR_SUP,
               "Particle field '%s' has no registered Eulerian scatter target.",
               particle_descriptor->canonical_name);
    PetscCheck(particle_descriptor->eulerian_scatter_target != FIELD_ID_INVALID,
               PETSC_COMM_SELF, PETSC_ERR_PLIB,
               "Particle field '%s' advertises scatter support without an Eulerian target.",
               particle_descriptor->canonical_name);
    particleFieldName = particle_descriptor->canonical_name;
    expected_dof = particle_descriptor->components;
    ierr = FieldGetDescriptor(particle_descriptor->eulerian_scatter_target, &eulerian_descriptor); CHKERRQ(ierr);
    PetscCheck(eulerian_descriptor->dof == expected_dof, PETSC_COMM_SELF, PETSC_ERR_PLIB,
               "Particle field '%s' has %d components but Eulerian target '%s' has %d.",
               particleFieldName, expected_dof, eulerian_descriptor->canonical_name, eulerian_descriptor->dof);
    ierr = FieldGetView(user, particle_descriptor->eulerian_scatter_target, &target_view); CHKERRQ(ierr);
    targetDM = target_view.dm;
 
    // --- 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, "%s", 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, "%s", 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, "%s", 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
                                                     particle_field_id,
                                                     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.

Definition at line 1888 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, ParticleFieldName(PARTICLE_FIELD_ID_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, ParticleFieldName(PARTICLE_FIELD_ID_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, PARTICLE_FIELD_ID_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");
    }
 
    // Additional scatterable particle fields require an explicit catalog entry
    // with a compatible persistent Eulerian target.
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished scattering specified particle fields.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ScatterParticleFieldToEulerField_Internal()

static PetscErrorCode ScatterParticleFieldToEulerField_Internal ( UserCtx *  user, ParticleFieldId  particle_field_id, DM  targetDM, PetscInt  expected_dof, Vec  eulerFieldAverageVec  )

static

Accumulate one particle field onto the Eulerian grid using the selected scatter stencil.

Definition at line 1715 of file interpolation.c.

{
    PetscErrorCode ierr;
    PetscInt       target_dof = 0;
    Vec            globalsumVec = NULL;
    Vec            localsumVec = NULL;
    char           msg[ERROR_MSG_BUFFER_SIZE]; // Buffer for formatted error messages
    const char    *particleFieldName = ParticleFieldName(particle_field_id);
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (!user || !user->swarm || !user->ParticleCount || !targetDM || !eulerFieldAverageVec)
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "NULL input provided to ScatterParticleFieldToEulerField_Internal.");
 
    ierr = DMDAGetInfo(targetDM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &target_dof, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    if (target_dof != expected_dof) {
        PetscSNPrintf(msg, sizeof(msg),
                      "Field '%s' expects DOF %d but targetDM reports DOF %d.",
                      particleFieldName, expected_dof, target_dof);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "%s", msg);
    }
    
    // --- Check if Particle Field Exists ---
    // Attempt a GetField call; if it fails, the field doesn't exist.
    // We let CHKERRQ handle the error directly if the field doesn't exist OR
    // we catch it specifically to provide a more tailored message.
 
    /*
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"Field %s being accessed to check existence \n",particleFieldName);
    ierr = DMSwarmGetField(user->swarm, particleFieldName, NULL, NULL, NULL);
    if (ierr) { // If GetField returns an error
        PetscSNPrintf(msg, sizeof(msg), "Particle field '%s' not found in DMSwarm for scattering.", particleFieldName);
        // Directly set the error, overwriting the one from GetField
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, msg);
    }
    ierr = DMSwarmRestoreField(user->swarm, particleFieldName, NULL, NULL, NULL);
    */
    
    // --- Setup Temporary Sum Vector ---
    ierr = VecDuplicate(eulerFieldAverageVec, &globalsumVec); CHKERRQ(ierr);
    ierr = VecSet(globalsumVec, 0.0); CHKERRQ(ierr);
    ierr = PetscSNPrintf(msg, sizeof(msg), "TempSum_%s", particleFieldName); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)globalsumVec, msg); CHKERRQ(ierr);
 
    // create local vector for accumulation
    ierr = DMGetLocalVector(targetDM, &localsumVec); CHKERRQ(ierr);
    ierr = VecSet(localsumVec, 0.0); CHKERRQ(ierr);  // Must be zeroed before accumulation
    ierr = PetscSNPrintf(msg, sizeof(msg), "LocalTempSum_%s", particleFieldName); CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)localsumVec, msg); CHKERRQ(ierr);
 
    // --- Accumulate ---
    // This will call DMSwarmGetField again. If it failed above, it will likely fail here too,
    // unless the error was cleared somehow between the check and here (unlikely).
    // If the check above was skipped (Option 1), this is where the error for non-existent
    // field will be caught by CHKERRQ.
    ierr = AccumulateParticleField(user->swarm, particle_field_id, targetDM, localsumVec); CHKERRQ(ierr);
 
    // --- Local to Global Sum ---
    ierr = DMLocalToGlobalBegin(targetDM, localsumVec, ADD_VALUES, globalsumVec); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(targetDM, localsumVec, ADD_VALUES, globalsumVec); CHKERRQ(ierr);
    // Return local vector to DM
    ierr = DMRestoreLocalVector(targetDM, &localsumVec); CHKERRQ(ierr);
 
    // Calculate the number of particles per cell.
    ierr = CalculateParticleCountPerCell(user); CHKERRQ(ierr);
    // --- Normalize ---
    ierr = NormalizeGridVectorByCount(user->da, user->ParticleCount, targetDM, globalsumVec, eulerFieldAverageVec); CHKERRQ(ierr);
 
    // --- Cleanup ---
    ierr = VecDestroy(&globalsumVec); CHKERRQ(ierr);
 
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: