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.
static PetscErrorCode	ScatterParticleFieldToEulerField_Internal (UserCtx *user, const char *particleFieldName, DM targetDM, PetscInt expected_dof, Vec eulerFieldAverageVec)
	Internal helper implementation: ScatterParticleFieldToEulerField_Internal().

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.

Definition at line 1486 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, "%s", 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.

Definition at line 1895 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, "%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
                                                     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.

Definition at line 1965 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:

Here is the caller graph for this function:

ScatterParticleFieldToEulerField_Internal()

static PetscErrorCode ScatterParticleFieldToEulerField_Internal ( UserCtx *  user, const char *  particleFieldName, DM  targetDM, PetscInt  expected_dof, Vec  eulerFieldAverageVec  )

static

Internal helper implementation: ScatterParticleFieldToEulerField_Internal().

Local to this translation unit.

Definition at line 1809 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
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (!user || !user->swarm || !user->ParticleCount || !particleFieldName || !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, particleFieldName, 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: