ParticleMotion.h File Reference

Header file for Particle Motion and migration related functions. More...

#include <petsc.h>
#include <petscdmswarm.h>
#include <stdbool.h>
#include <petscsys.h>
#include <math.h>
#include "variables.h"
#include "logging.h"
#include "walkingsearch.h"

Include dependency graph for ParticleMotion.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	GenerateGaussianNoise (PetscRandom rnd, PetscReal *n1, PetscReal *n2)
	Generates two independent standard normal random variables N(0,1) using the Box-Muller transform.
PetscErrorCode	CalculateBrownianDisplacement (UserCtx *user, PetscReal diff_eff, Cmpnts *displacement)
	Calculates the stochastic displacement vector (Brownian motion) for a single particle.
PetscErrorCode	UpdateParticlePosition (UserCtx *user, Particle *particle)
	Updates a particle's position based on its velocity and the timestep dt (stored in user->dt).
PetscErrorCode	UpdateAllParticlePositions (UserCtx *user)
	Loops over all local particles in the DMSwarm, updating their positions based on velocity and the global timestep user->dt.
PetscErrorCode	CheckAndRemoveOutOfBoundsParticles (UserCtx *user, PetscInt *removedCountLocal, PetscInt *removedCountGlobal, const BoundingBox *bboxlist)
	Checks for particles outside the physical domain boundaries and removes them using DMSwarmRemovePointAtIndex.
PetscErrorCode	CheckAndRemoveLostParticles (UserCtx *user, PetscInt *removedCountLocal, PetscInt *removedCountGlobal, PetscReal *removedScalarSumGlobal)
	Removes particles that have been definitively flagged as LOST by the location algorithm.
PetscErrorCode	DefineBasicMigrationPattern (UserCtx *user)
	Defines the basic migration pattern for particles within the swarm.
PetscErrorCode	PerformBasicMigration (UserCtx *user)
	Performs the basic migration of particles based on the defined migration pattern.
PetscErrorCode	IdentifyMigratingParticles (UserCtx *user, const BoundingBox *bboxlist, MigrationInfo **migrationList, PetscInt *migrationCount, PetscInt *listCapacity)
	Identifies particles leaving the local bounding box and finds their target neighbor rank.
PetscErrorCode	SetMigrationRanks (UserCtx *user, const MigrationInfo *migrationList, PetscInt migrationCount)
	Writes migration destinations into the DMSwarm rank field for marked particles.
PetscErrorCode	PerformMigration (UserCtx *user)
	Performs particle migration based on the pre-populated DMSwarmPICField_rank field.
PetscErrorCode	CalculateParticleCountPerCell (UserCtx *user)
	Counts particles in each cell of the DMDA 'da' and stores the result in user->ParticleCount.
PetscErrorCode	ResizeSwarmGlobally (DM swarm, PetscInt N_target)
	Resizes a swarm collectively to a target global particle count.
PetscErrorCode	PreCheckAndResizeSwarm (UserCtx *user, PetscInt ti, const char *ext)
	Checks particle count in the reference file and resizes the swarm if needed.
PetscErrorCode	PerformSingleParticleMigrationCycle (UserCtx *user, const BoundingBox *bboxlist, MigrationInfo **migrationList_p, PetscInt *migrationCount_p, PetscInt *migrationListCapacity_p, PetscReal currentTime, PetscInt step, const char *migrationCycleName, PetscInt *globalMigrationCount_out)
	Performs one full cycle of particle migration: identify, set ranks, and migrate.
PetscErrorCode	ReinitializeParticlesOnInletSurface (UserCtx *user, PetscReal currentTime, PetscInt step)
	Re-initializes the positions of particles currently on this rank if this rank owns part of the designated inlet surface.
PetscErrorCode	GetLocalPIDSnapshot (const PetscInt64 pid_field[], PetscInt n_local, PetscInt64 **pids_snapshot_out)
	Creates a sorted snapshot of all Particle IDs (PIDs) from a raw data array.
PetscErrorCode	AddToMigrationList (MigrationInfo **migration_list_p, PetscInt *capacity_p, PetscInt *count_p, PetscInt particle_local_idx, PetscMPIInt destination_rank)
	Safely adds a new migration task to a dynamically sized list.
PetscErrorCode	FlagNewcomersForLocation (DM swarm, PetscInt n_local_before, const PetscInt64 pids_before[])
	Identifies newly arrived particles after migration and flags them for a location search.
PetscErrorCode	MigrateRestartParticlesUsingCellID (UserCtx *user)
	Fast-path migration for restart particles using preloaded Cell IDs.
PetscErrorCode	LocateAllParticlesInGrid (UserCtx *user, BoundingBox *bboxlist)
	Orchestrates the complete particle location and migration process for one timestep.
PetscErrorCode	ResetAllParticleStatuses (UserCtx *user)
	Marks all local particles as NEEDS_LOCATION for the next settlement pass.

Detailed Description

Header file for Particle Motion and migration related functions.

This file contains declarations of functions responsible for moving and migrating particle swarms within a simulation using PETSc's DMSwarm.

Definition in file ParticleMotion.h.

Function Documentation

GenerateGaussianNoise()

PetscErrorCode GenerateGaussianNoise	(	PetscRandom	rnd,
		PetscReal *	n1,
		PetscReal *	n2
	)

Generates two independent standard normal random variables N(0,1) using the Box-Muller transform.

Parameters

[in]	rnd	The PETSc Random context (Uniform [0,1)).
[out]	n1	First Gaussian number.
[out]	n2	Second Gaussian number.

Returns

PetscErrorCode

Generates two independent standard normal random variables N(0,1) using the Box-Muller transform.

Local to this translation unit.

Definition at line 17 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscScalar    val1, val2;
    PetscReal      u1, u2;
    PetscReal      magnitude, theta;
 
    PetscFunctionBeginUser;
 
    // 1. Get two independent uniform random numbers from the generator
    // PetscRandomGetValue returns a PetscScalar (which might be complex).
    // We take the Real part to ensure this works in both Real and Complex builds.
    ierr = PetscRandomGetValue(rnd, &val1); CHKERRQ(ierr);
    ierr = PetscRandomGetValue(rnd, &val2); CHKERRQ(ierr);
 
    u1 = PetscRealPart(val1);
    u2 = PetscRealPart(val2);
 
    // 2. Safety Check: log(0) is undefined (infinity).
    // If the RNG returns exactly 0.0, bump it to a tiny epsilon.
    if (u1 <= 0.0) u1 = 1.0e-14;
 
    // 3. Box-Muller Transform
    // Formula: R = sqrt(-2 * ln(u1)), Theta = 2 * PI * u2
    magnitude = PetscSqrtReal(-2.0 * PetscLogReal(u1));
    theta     = 2.0 * PETSC_PI * u2; 
 
    // 4. Calculate independent Normal variables
    *n1 = magnitude * PetscCosReal(theta);
    *n2 = magnitude * PetscSinReal(theta);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

CalculateBrownianDisplacement()

PetscErrorCode CalculateBrownianDisplacement	(	UserCtx *	user,
		PetscReal	diff_eff,
		Cmpnts *	displacement
	)

Calculates the stochastic displacement vector (Brownian motion) for a single particle.

Equation: dX_stoch = sqrt(2 * Gamma_eff * dt) * N(0,1)

Parameters

[in]	user	Pointer to UserCtx (access to dt and BrownianMotionRNG).
[in]	diff_eff	The effective diffusivity (Gamma + Gamma_t) at the particle's location.
[out]	displacement	Pointer to a Cmpnts struct to store the resulting (dx, dy, dz).

Returns

PetscErrorCode

Calculates the stochastic displacement vector (Brownian motion) for a single particle.

Local to this translation unit.

Definition at line 57 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscReal      dt = user->simCtx->dt;
    PetscReal      sigma;
    PetscReal      n_x, n_y, n_z, gaussian_dummy;
 
    PetscFunctionBeginUser;
 
    // 1. Initialize output to zero for safety
    displacement->x = 0.0;
    displacement->y = 0.0;
    displacement->z = 0.0;
 
    // 2. Physical check: Diffusivity cannot be negative. 
    // If 0, there is no Brownian motion.
    if (diff_eff <= 1.0e-12) {
        PetscFunctionReturn(0);
    }
 
    // 3. Calculate the Scaling Factor (Standard Deviation)
    // Formula: sigma = sqrt(2 * D * dt)
    // Note: dt is inside the root because variance scales linearly with time.
    sigma = PetscSqrtReal(2.0 * diff_eff * dt);
 
    // 4. Generate 3 Independent Gaussian Random Numbers
    // GenerateGaussianNoise produces 2 numbers at a time. We call it twice.
    
    // Get noise for X and Y
    ierr = GenerateGaussianNoise(user->simCtx->BrownianMotionRNG, &n_x, &n_y); CHKERRQ(ierr);
    
    // Get noise for Z (second sample is intentionally discarded here).
    ierr = GenerateGaussianNoise(user->simCtx->BrownianMotionRNG, &n_z, &gaussian_dummy); CHKERRQ(ierr);
 
    // 5. Calculate final stochastic displacement
    displacement->x = sigma * n_x;
    displacement->y = sigma * n_y;
    displacement->z = sigma * n_z;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

UpdateParticlePosition()

PetscErrorCode UpdateParticlePosition	(	UserCtx *	user,
		Particle *	particle
	)

Updates a particle's position based on its velocity and the timestep dt (stored in user->dt).

Parameters

[in]	user	Pointer to your UserCtx (must contain user->dt).
[in,out]	particle	Pointer to the particle struct (contains, pos,vel,diffusivity etc).

Returns

PetscErrorCode Returns 0 on success, or an error code on failure.

Updates a particle's position based on its velocity and the timestep dt (stored in user->dt).

Local to this translation unit.

Definition at line 105 of file ParticleMotion.c.

{
  PetscFunctionBeginUser; // PETSc macro for error/stack tracing
  PROFILE_FUNCTION_BEGIN;
 
  PetscErrorCode ierr;
  PetscReal  dt = user->simCtx->dt;
  Cmpnts brownian_disp;
 
  // 2. Calculate the stochastic kick
  ierr = CalculateBrownianDisplacement(user,particle->diffusivity, &brownian_disp); CHKERRQ(ierr);
 
  // --- Update Position ---
  // X_new = X_old + ((U_convection + U_diffusivitygradient) * dt) + dX_brownian
  
  particle->loc.x += ((particle->vel.x + particle->diffusivitygradient.x) * dt) + brownian_disp.x;
  particle->loc.y += ((particle->vel.y + particle->diffusivitygradient.y) * dt) + brownian_disp.y;
  particle->loc.z += ((particle->vel.z + particle->diffusivitygradient.z) * dt) + brownian_disp.z;
 
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

UpdateAllParticlePositions()

PetscErrorCode UpdateAllParticlePositions

(

UserCtx *

user

)

Loops over all local particles in the DMSwarm, updating their positions based on velocity and the global timestep user->dt.

Parameters

[in,out]

user

Pointer to UserCtx (must contain dt).

Returns

PetscErrorCode Returns 0 on success, or an error code on failure.

Loops over all local particles in the DMSwarm, updating their positions based on velocity and the global timestep user->dt.

Local to this translation unit.

Definition at line 134 of file ParticleMotion.c.

{
  PetscErrorCode ierr;
  DM swarm = user->swarm;
  PetscInt       nLocal, p;
  PetscReal        *pos = NULL;
  PetscReal        *vel = NULL;
  PetscReal        *diffusivity = NULL;
  Cmpnts           *diffusivitygradient = NULL; 
  PetscReal        *psi = NULL;
  PetscReal        *weights = NULL;
  PetscInt         *cell = NULL;
  PetscInt         *status = NULL;
  PetscInt64       *pid = NULL;
  PetscMPIInt rank;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 
  PetscFunctionBeginUser;  // PETSc macro for error/stack tracing
  
  PROFILE_FUNCTION_BEGIN;
 
  // 1) Get the number of local particles
  ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
  if (nLocal == 0) {
    LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] No particles to move/transport. \n",rank);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);   // nothing to do, no fields held 
  }
  // 2) Access the "position" and "velocity" fields
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), NULL, NULL, (void**)&vel); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_DIFFUSIVITY), NULL, NULL, (void**)&diffusivity); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT), NULL, NULL, (void**)&diffusivitygradient); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI), NULL, NULL, (void**)&psi); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT), NULL, NULL, (void**)&weights); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cell); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status); CHKERRQ(ierr);
  ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void**)&pid); CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG," [Rank %d] No.of Particles to update: %" PetscInt_FMT ".\n",rank,nLocal);
 
  // 3) Loop over all local particles, updating each position by velocity * dt
  for (p = 0; p < nLocal; p++) {
    // update temporary particle struct
    Particle particle;
 
    // Unpack: Use the helper to read from swarm arrays into the particle struct
    ierr = UnpackSwarmFields(p, pid, weights, pos, cell, vel, status, diffusivity, diffusivitygradient, psi, &particle); CHKERRQ(ierr);
 
    // Update position based on velocity and Brownian motion
    ierr = UpdateParticlePosition(user, &particle); CHKERRQ(ierr);
    
    // Update swarm fields
    ierr = UpdateSwarmFields(p, &particle, pos, vel, weights, cell, status, diffusivity, diffusivitygradient, psi); CHKERRQ(ierr);
  }
 
  // 4) Restore the fields
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY), NULL, NULL, (void**)&vel); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_DIFFUSIVITY), NULL, NULL, (void**)&diffusivity); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_DIFFUSIVITY_GRADIENT), NULL, NULL, (void**)&diffusivitygradient); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI), NULL, NULL, (void**)&psi); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT), NULL, NULL, (void**)&weights); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cell); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status); CHKERRQ(ierr);
  ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void**)&pid); CHKERRQ(ierr);
 
 
  LOG_ALLOW(LOCAL,LOG_DEBUG,"Particle moved/transported successfully on Rank %d.\n",rank);
  
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CheckAndRemoveOutOfBoundsParticles()

PetscErrorCode CheckAndRemoveOutOfBoundsParticles	(	UserCtx *	user,
		PetscInt *	removedCountLocal,
		PetscInt *	removedCountGlobal,
		const BoundingBox *	bboxlist
	)

Checks for particles outside the physical domain boundaries and removes them using DMSwarmRemovePointAtIndex.

This function iterates through all particles local to the current MPI rank. It checks if a particle's position (x, y, or z) is outside the specified physical domain boundaries [xMin, xMax], [yMin, yMax], [zMin, zMax].

If a particle is found out of bounds, it is removed using DMSwarmRemovePointAtIndex. NOTE: Removing points changes the indices of subsequent points in the iteration. Therefore, it's crucial to iterate BACKWARDS or carefully manage indices after a removal. Iterating backwards is generally safer.

Parameters

	user	Pointer to the UserCtx structure.
[out]	removedCountLocal	Pointer to store the number of particles removed on this rank.
[out]	removedCountGlobal	Pointer to store the total number of particles removed across all ranks.
[in]	bboxlist	An array of BoundingBox structures for ALL MPI ranks, indexed 0 to (size-1). This array must be up-to-date and available on all ranks.

Returns

PetscErrorCode 0 on success, non-zero on failure.

Checks for particles outside the physical domain boundaries and removes them using DMSwarmRemovePointAtIndex.

Local to this translation unit.

Definition at line 228 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    PetscInt       nLocalInitial;
    PetscReal      *pos_p = NULL;
    PetscInt64     *pid_p = NULL; // For better logging
    PetscInt       local_removed_count = 0;
    PetscMPIInt    global_removed_count_mpi = 0;
    PetscMPIInt    rank, size;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_INFO, "[Rank %d] Checking for out-of-bounds particles...", rank);
 
    // Initialize output parameters to ensure clean state
    *removedCountLocal = 0;
    if (removedCountGlobal) *removedCountGlobal = 0;
 
    ierr = DMSwarmGetLocalSize(swarm, &nLocalInitial); CHKERRQ(ierr);
 
    // Only proceed if there are particles to check on this rank.
    // All ranks will still participate in the final collective MPI_Allreduce.
    if (nLocalInitial > 0) {
        // Get access to swarm fields once before the loop begins.
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),    NULL, NULL, (void **)&pos_p); CHKERRQ(ierr);
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void **)&pid_p); CHKERRQ(ierr);
 
        // --- Iterate BACKWARDS to handle index changes safely during removal ---
        for (PetscInt p = nLocalInitial - 1; p >= 0; p--) {
            PetscBool isInsideAnyBox = PETSC_FALSE;
            Cmpnts current_pos = {pos_p[3*p + 0], pos_p[3*p + 1], pos_p[3*p + 2]};
 
            // Check if the particle is inside ANY of the rank bounding boxes
            for (PetscMPIInt proc = 0; proc < size; proc++) {
                if (IsParticleInBox(&bboxlist[proc], &current_pos)) {
                    isInsideAnyBox = PETSC_TRUE;
                    break; // Particle is inside a valid domain, stop checking.
                }
            }
 
            if (!isInsideAnyBox) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Removing out-of-bounds particle [PID %lld] at local index %d. Pos: (%g, %g, %g)\n",
                          rank, (long long)pid_p[p], p, current_pos.x, current_pos.y, current_pos.z);
 
                // --- Safe Removal Pattern: Restore -> Remove -> Reacquire ---
                // This is the fix for the double-restore bug. Pointers are managed carefully
                // within this block and then restored cleanly after the loop.
 
                // 1. Restore all fields BEFORE modifying the swarm structure. This invalidates pos_p and pid_p.
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),    NULL, NULL, (void **)&pos_p); CHKERRQ(ierr);
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void **)&pid_p); CHKERRQ(ierr);
 
                // 2. Remove the particle at the current local index 'p'.
                ierr = DMSwarmRemovePointAtIndex(swarm, p); CHKERRQ(ierr);
                local_removed_count++;
 
                // 3. After removal, re-acquire pointers ONLY if the loop is not finished.
                PetscInt nLocalCurrent;
                ierr = DMSwarmGetLocalSize(swarm, &nLocalCurrent); CHKERRQ(ierr);
 
                if (nLocalCurrent > 0 && p > 0) { // Check if there are particles left AND iterations left
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),    NULL, NULL, (void **)&pos_p); CHKERRQ(ierr);
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void **)&pid_p); CHKERRQ(ierr);
                } else {
                    // All remaining particles were removed OR this was the last particle (p=0).
                    // Invalidate pointers to prevent the final restore call and exit the loop.
                    pos_p = NULL;
                    pid_p = NULL;
                    break;
                }
            }
        } // End of backwards loop
 
        // At the end, restore any valid pointers. This handles three cases:
        // 1. No particles were removed: restores the original pointers.
        // 2. Particles were removed mid-loop: restores the pointers from the last re-acquisition.
        // 3. All particles were removed: pointers are NULL, so nothing is done.
        if (pos_p) { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),    NULL, NULL, (void **)&pos_p); CHKERRQ(ierr); }
        if (pid_p) { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID), NULL, NULL, (void **)&pid_p); CHKERRQ(ierr); }
    } // End of if (nLocalInitial > 0)
 
    PetscInt nLocalFinal;
    ierr = DMSwarmGetLocalSize(swarm, &nLocalFinal); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_INFO, "[Rank %d] Finished removing %d out-of-bounds particles. Final local size: %d.\n", rank, local_removed_count, nLocalFinal);
 
    // --- Synchronize counts across all ranks ---
    *removedCountLocal = local_removed_count;
    if (removedCountGlobal) {
        ierr = MPI_Allreduce(&local_removed_count, &global_removed_count_mpi, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)swarm)); CHKERRQ(ierr);
        *removedCountGlobal = global_removed_count_mpi;
        // Use a synchronized log message so only one rank prints the global total.
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Removed %d out-of-bounds particles globally.\n", rank, *removedCountGlobal);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CheckAndRemoveLostParticles()

PetscErrorCode CheckAndRemoveLostParticles	(	UserCtx *	user,
		PetscInt *	removedCountLocal,
		PetscInt *	removedCountGlobal,
		PetscReal *	removedScalarSumGlobal
	)

Removes particles that have been definitively flagged as LOST by the location algorithm.

This function is the designated cleanup utility. It should be called after the LocateAllParticlesInGrid orchestrator has run and every particle's status has been definitively determined.

It iterates through all locally owned particles and checks their DMSwarm_location_status field. If a particle's status is LOST, it is permanently removed from the simulation using DMSwarmRemovePointAtIndex.

This approach centralizes the removal logic, making the DMSwarm_location_status the single source of truth for a particle's validity, which is more robust than relying on secondary geometric checks (like bounding boxes).

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the swarm.
[out]	removedCountLocal	Pointer to store the number of particles removed on this rank.
[out]	removedCountGlobal	Pointer to store the total number of particles removed across all ranks.
[out]	removedScalarSumGlobal	Optional, may be NULL: the sum of Psi over every particle removed on any rank. With a 0/1 label and micromixing off, it counts the removed particles that carried the label.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Removes particles that have been definitively flagged as LOST by the location algorithm.

Local to this translation unit.

Definition at line 336 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    PetscInt       nLocalInitial;
    PetscInt       *status_p = NULL;
    PetscInt64     *pid_p = NULL; // For better logging
    PetscReal      *pos_p = NULL; // For better logging
    PetscReal      *psi_p = NULL;
    PetscReal      local_removed_psi = 0.0;
    PetscInt       local_removed_count = 0;
    PetscMPIInt    global_removed_count_mpi = 0;
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Checking for and removing LOST particles...\n", rank);
 
    // Initialize output parameters to ensure clean state
    *removedCountLocal = 0;
    if (removedCountGlobal) *removedCountGlobal = 0;
    if (removedScalarSumGlobal) *removedScalarSumGlobal = 0.0;
 
    ierr = DMSwarmGetLocalSize(swarm, &nLocalInitial); CHKERRQ(ierr);
 
    // Only proceed if there are particles to check on this rank.
    // All ranks will still participate in the final collective MPI_Allreduce.
    if (nLocalInitial > 0) {
        // Get access to all swarm fields once before the loop begins.
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void **)&status_p); CHKERRQ(ierr);
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void **)&pid_p);    CHKERRQ(ierr);
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void **)&pos_p);    CHKERRQ(ierr);
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI),                     NULL, NULL, (void **)&psi_p);    CHKERRQ(ierr);
 
        // --- Iterate BACKWARDS to handle index changes safely during removal ---
        for (PetscInt p = nLocalInitial - 1; p >= 0; p--) {
            if (status_p[p] == LOST) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Removing LOST particle [PID %lld] at local index %d. Position: (%.4f, %.4f, %.4f).\n",
                          rank, (long long)pid_p[p], p, pos_p[3*p], pos_p[3*p+1], pos_p[3*p+2]);
 
                // --- Safe Removal Pattern: Restore -> Remove -> Reacquire ---
                // This is the fix for the double-restore bug. Pointers are managed carefully
                // within this block and then restored cleanly after the loop.
 
                // The scalar a removed particle carries leaves the domain with it.
                local_removed_psi += psi_p[p];
 
                // 1. Restore all fields BEFORE modifying the swarm structure. This invalidates all pointers.
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI),             NULL, NULL, (void **)&psi_p);    CHKERRQ(ierr);
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void **)&status_p); CHKERRQ(ierr);
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void **)&pid_p);    CHKERRQ(ierr);
                ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void **)&pos_p);    CHKERRQ(ierr);
 
                // 2. Remove the particle at the current local index 'p'.
                ierr = DMSwarmRemovePointAtIndex(swarm, p); CHKERRQ(ierr);
                local_removed_count++;
 
                // 3. After removal, re-acquire pointers ONLY if the loop is not finished.
                PetscInt nLocalCurrent;
                ierr = DMSwarmGetLocalSize(swarm, &nLocalCurrent); CHKERRQ(ierr);
 
                if (nLocalCurrent > 0 && p > 0) { // Check if there are particles left AND iterations left
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void **)&status_p); CHKERRQ(ierr);
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void **)&pid_p);    CHKERRQ(ierr);
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void **)&pos_p);    CHKERRQ(ierr);
                    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI),                     NULL, NULL, (void **)&psi_p);    CHKERRQ(ierr);
                } else {
                    // All remaining particles were removed OR this was the last particle (p=0).
                    // Invalidate pointers to prevent the final restore call and exit the loop.
                    status_p = NULL;
                    pid_p = NULL;
                    pos_p = NULL;
                    psi_p = NULL;
                    break;
                }
            }
        } // End of backwards loop
 
        // At the end, restore any valid pointers. This handles three cases:
        // 1. No particles were removed: restores the original pointers.
        // 2. Particles were removed mid-loop: restores the pointers from the last re-acquisition.
        // 3. All particles were removed: pointers are NULL, so nothing is done.
        if (status_p) { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void **)&status_p); CHKERRQ(ierr); }
        if (pid_p)    { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void **)&pid_p);    CHKERRQ(ierr); }
        if (pos_p)    { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void **)&pos_p);    CHKERRQ(ierr); }
        if (psi_p)    { ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PSI),                     NULL, NULL, (void **)&psi_p);    CHKERRQ(ierr); }
    } // End of if (nLocalInitial > 0)
 
    PetscInt nLocalFinal;
    ierr = DMSwarmGetLocalSize(swarm, &nLocalFinal); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d: Finished removing %d LOST particles. Final local size: %d.\n", rank, local_removed_count, nLocalFinal);
 
    // --- Synchronize counts across all ranks ---
    *removedCountLocal = local_removed_count;
    if (removedCountGlobal) {
        ierr = MPI_Allreduce(&local_removed_count, &global_removed_count_mpi, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)swarm)); CHKERRQ(ierr);
        *removedCountGlobal = global_removed_count_mpi;
        // Use a synchronized log message so only one rank prints the global total.
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Removed %d LOST particles globally.\n", rank, *removedCountGlobal);
    }
    if (removedScalarSumGlobal) {
        ierr = MPI_Allreduce(&local_removed_psi, removedScalarSumGlobal, 1, MPIU_REAL, MPI_SUM,
                             PetscObjectComm((PetscObject)swarm)); CHKERRMPI(ierr);
    }
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DefineBasicMigrationPattern()

PetscErrorCode DefineBasicMigrationPattern

(

UserCtx *

user

)

Defines the basic migration pattern for particles within the swarm.

This function establishes the migration pattern that dictates how particles move between different MPI ranks in the simulation. It initializes a migration list where each particle is assigned a target rank based on predefined conditions. The migration pattern can be customized to implement various migration behaviors.

Parameters

[in,out]

user

Pointer to the UserCtx structure containing simulation context.

Returns

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

PerformBasicMigration()

PetscErrorCode PerformBasicMigration

(

UserCtx *

user

)

Performs the basic migration of particles based on the defined migration pattern.

This function updates the positions of particles within the swarm by migrating them to target MPI ranks as specified in the migration list. It handles the migration process by setting the 'DMSwarm_rank' field for each particle and invokes the DMSwarm migration mechanism to relocate particles across MPI processes. After migration, it cleans up allocated resources and ensures synchronization across all MPI ranks.

Parameters

[in,out]

user

Pointer to the UserCtx structure containing simulation context.

Returns

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

IdentifyMigratingParticles()

PetscErrorCode IdentifyMigratingParticles	(	UserCtx *	user,
		const BoundingBox *	bboxlist,
		MigrationInfo **	migrationList,
		PetscInt *	migrationCount,
		PetscInt *	listCapacity
	)

Identifies particles leaving the local bounding box and finds their target neighbor rank.

Iterates local particles, checks against local bounding box. If outside, checks the pre-computed immediate neighbors (user->neighbors) using the global bboxlist to see if the particle landed in one of them. Populates the migrationList. Does NOT handle particles leaving the global domain (assumes CheckAndRemove was called).

Parameters

user	Pointer to the UserCtx (contains local bbox and neighbors).
bboxlist	Array of BoundingBox structs for all ranks (for checking neighbor boxes).
migrationList	Pointer to an array of MigrationInfo structs (output, allocated/reallocated by this func).
migrationCount	Pointer to the number of particles marked for migration (output).
listCapacity	Pointer to the current allocated capacity of migrationList (in/out).

Returns

PetscErrorCode 0 on success, non-zero on failure.

SetMigrationRanks()

PetscErrorCode SetMigrationRanks	(	UserCtx *	user,
		const MigrationInfo *	migrationList,
		PetscInt	migrationCount
	)

Writes migration destinations into the DMSwarm rank field for marked particles.

This helper consumes the migration list produced by IdentifyMigratingParticles and updates each selected particle's destination rank so that a subsequent PerformMigration call can transfer ownership correctly.

Parameters

[in,out]	user	Context containing the swarm and migration rank field.
[in]	migrationList	Array of migration directives (local index + destination rank).
[in]	migrationCount	Number of valid entries in migrationList.

Returns

PetscErrorCode 0 on success.

Writes migration destinations into the DMSwarm rank field for marked particles.

Local to this translation unit.

Definition at line 456 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    PetscInt       p_idx;
    PetscInt      *rankField = NULL; // Field storing target rank
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // Ensure the migration rank field exists
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_RANK), NULL, NULL, (void **)&rankField); CHKERRQ(ierr);
 
    // Set the target rank for migrating particles
    for(p_idx = 0; p_idx < migrationCount; ++p_idx) {
        rankField[migrationList[p_idx].local_index] = migrationList[p_idx].target_rank;
    }
 
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_RANK), NULL, NULL, (void **)&rankField); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PerformMigration()

PetscErrorCode PerformMigration

(

UserCtx *

user

)

Performs particle migration based on the pre-populated DMSwarmPICField_rank field.

Assumes SetMigrationRanks has already been called to mark particles with their target ranks. Calls DMSwarmMigrate to execute the communication and removal of un-migrated particles.

Parameters

user	Pointer to the UserCtx structure containing the swarm.

Returns

PetscErrorCode 0 on success, non-zero on failure.

Performs particle migration based on the pre-populated DMSwarmPICField_rank field.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/ParticleMotion.h.

See also

PerformMigration()

Definition at line 489 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    DM swarm = user->swarm;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting DMSwarmMigrate...\n", rank);
 
    // Perform the migration - PETSC_TRUE removes particles that fail to land
    // in a valid cell on the target rank (or were marked with an invalid rank).
    ierr = DMSwarmMigrate(swarm, PETSC_TRUE); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Migration complete.\n", rank);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller 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.

Assumes user->ParticleCount is a pre-allocated global vector associated with user->da and initialized to zero before calling this function (though it resets it internally). Assumes particle 'DMSwarm_CellID' field contains local cell indices.

Parameters

[in,out]

user

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

Returns

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

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

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/logging.h.

See also

CalculateParticleCountPerCell()

Definition at line 521 of file ParticleMotion.c.

                                                            {
    PetscErrorCode ierr;
    DM             da = user->da;
    DM             swarm = user->swarm;
    Vec            countVec = user->ParticleCount;
    Vec            localcountVec = user->lParticleCount;
    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 %" PetscInt_FMT ".", count_dof);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", msg);
    }
 
    // --- Zero the local count vector ---
    ierr = VecSet(localcountVec, 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,ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void **)&global_cell_id_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm,ParticleFieldName(PARTICLE_FIELD_ID_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, localcountVec, &count_arr_3d); CHKERRQ(ierr);
 
    // Get local owned + ghosted range for writing into ghost slots.
    PetscInt gxs, gys, gzs, gxm, gym, gzm;
    ierr = DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gxm, &gym, &gzm); CHKERRQ(ierr);
 
    // --- Accumulate Counts Locally ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "CalculateParticleCountPerCell (Rank %d): Processing %" PetscInt_FMT " local particles using GLOBAL CellIDs.\n",rank,nlocal);
    for (p = 0; p < nlocal; p++) {
        // Read the GLOBAL indices stored for this particle
        PetscInt i_geom = global_cell_id_arr[p * 3 + 0]; // Global i index
        PetscInt j_geom = global_cell_id_arr[p * 3 + 1]; // Global j index
        PetscInt 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=%" PetscInt_FMT ", PID = %" PetscInt64_FMT ": (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n",
                       rank, p, PID_arr[p], i, j, k);
    
        // Check if the global index (i,j,k) falls within the local + ghost range
        if (i >= gxs && i < gxs + gxm  && 
            j >= gys && j < gys + gym  && // Adjust based on actual ghost width
            k >= gzs && k < gzs + gzm  )   // 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)
              // note: this is not necessarily an error if the particle is legitimately outside the local+ghost region
              LOG_ALLOW(LOCAL, LOG_VERBOSE,
                        "(Rank %d): Skipping particle %" PetscInt64_FMT " with global CellID (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") - likely outside local+ghost range.\n",
                        rank, PID_arr[p], i, j, k);
    }
    }
    LOG_ALLOW(LOCAL, LOG_DEBUG, "(Rank %d): Local counting finished. Processed %" PetscInt_FMT " particles locally.\n", rank, particles_counted_locally);
 
    // --- Restore Access ---
    ierr = DMDAVecRestoreArray(da, localcountVec, &count_arr_3d); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm,ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void **)&global_cell_id_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm,ParticleFieldName(PARTICLE_FIELD_ID_PID),NULL,NULL,(void **)&PID_arr);CHKERRQ(ierr);
 
    // --- Assemble Global Vector ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling global ParticleCount vector.\n");
    ierr = VecZeroEntries(countVec); CHKERRQ(ierr); // Ensure global vector is zeroed before accumulation
    ierr = DMLocalToGlobalBegin(da, localcountVec, ADD_VALUES, countVec); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(da, localcountVec, ADD_VALUES, countVec); CHKERRQ(ierr);
    /* 
     * OPTIONAL: Synchronize Ghosts for Stencil Operations
     * If a future function needs to read ParticleCount from neighbor cells (e.g., density smoothing 
     * or gradient calculations), uncomment the following lines to update the ghost slots 
     * in user->lParticleCount with the final summed values.
     * 
     ierr = UpdateLocalGhosts(user, FIELD_ID_PARTICLE_COUNT); 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 call graph for this function:

Here is the caller graph for this function:

ResizeSwarmGlobally()

PetscErrorCode ResizeSwarmGlobally	(	DM	swarm,
		PetscInt	N_target
	)

Resizes a swarm collectively to a target global particle count.

The target is divided by quotient and remainder so every rank receives either floor(N_target / nranks) or one additional entry. Resizing establishes the storage layout; callers initialize or overwrite particle fields afterwards.

Parameters

[in,out]	swarm	Swarm object to resize.
[in]	N_target	Target global particle count.

Returns

PetscErrorCode 0 on success.

Resizes a swarm collectively to a target global particle count.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/ParticleMotion.h.

See also

ResizeSwarmGlobally()

Definition at line 690 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscInt       N_current, N_final, nlocal_current, nlocal_target;
    PetscMPIInt    rank, size;
    MPI_Comm       comm;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = PetscObjectGetComm((PetscObject)swarm, &comm); CHKERRQ(ierr);
    ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr);
    PetscCheck(N_target >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE,
               "Target swarm size must be nonnegative; got %" PetscInt_FMT ".", N_target);
    ierr = DMSwarmGetSize(swarm, &N_current); CHKERRQ(ierr);
    ierr = DMSwarmGetLocalSize(swarm, &nlocal_current); CHKERRQ(ierr);
    nlocal_target = N_target / size + (rank < N_target % size ? 1 : 0);
 
    if (nlocal_current != nlocal_target) {
        LOG_ALLOW(LOCAL, LOG_DEBUG,
                  "Rank %d: resizing local swarm share from %" PetscInt_FMT
                  " to %" PetscInt_FMT ".\n",
                  rank, nlocal_current, nlocal_target);
        ierr = DMSwarmSetLocalSizes(swarm, nlocal_target, -1); CHKERRQ(ierr);
    }
 
    // Verify final size
    ierr = DMSwarmGetSize(swarm, &N_final); CHKERRQ(ierr);
    if (N_final != N_target) {
        SETERRQ(comm, PETSC_ERR_PLIB,
                "Failed to resize swarm: expected %" PetscInt_FMT
                " particles, got %" PetscInt_FMT, N_target, N_final);
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Swarm resized from %" PetscInt_FMT " to %" PetscInt_FMT " particles.\n",
              N_current, N_final);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PreCheckAndResizeSwarm()

PetscErrorCode PreCheckAndResizeSwarm	(	UserCtx *	user,
		PetscInt	ti,
		const char *	ext
	)

Checks particle count in the reference file and resizes the swarm if needed.

Reads the specified field file (e.g., position) into a temporary Vec to determine the number of particles (N_file) represented in that file for the given timestep. Compares N_file with the current swarm size (N_current). If they differ, resizes the swarm globally (adds or removes particles) to match N_file. The resized population is balanced across the communicator before field input.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the DMSwarm.
[in]	ti	Time index for constructing the file name.
[in]	ext	File extension (e.g., "dat").

Returns

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

Checks particle count in the reference file and resizes the swarm if needed.

Local to this translation unit.

Definition at line 736 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscInt       N_file = 0;
    PetscInt       N_current = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    (void)ext;
    ierr = ReadCheckpointParticleCount(user, ti, &N_file); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Committed checkpoint step %d records %d particles.\n", ti, N_file);
 
 
    // --- Now all ranks have the correct N_file, compare and resize if needed ---
    ierr = DMSwarmGetSize(user->swarm, &N_current); CHKERRQ(ierr);
 
    if (N_file != N_current) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Swarm size %d differs from file size %d. Resizing swarm globally.\n", N_current, N_file);
        ierr = ResizeSwarmGlobally(user->swarm, N_file); CHKERRQ(ierr);
    } else {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Swarm size (%d) already matches file size. No resize needed.\n", N_current);
    }
    
    // Also update the context
    user->simCtx->np = N_file;
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PerformSingleParticleMigrationCycle()

PetscErrorCode PerformSingleParticleMigrationCycle	(	UserCtx *	user,
		const BoundingBox *	bboxlist,
		MigrationInfo **	migrationList_p,
		PetscInt *	migrationCount_p,
		PetscInt *	migrationListCapacity_p,
		PetscReal	currentTime,
		PetscInt	step,
		const char *	migrationCycleName,
		PetscInt *	globalMigrationCount_out
	)

Performs one full cycle of particle migration: identify, set ranks, and migrate.

This function encapsulates the three main steps of migrating particles between MPI ranks:

1. Identify particles on the local rank that need to move based on their current positions and the domain decomposition (bboxlist).
2. Determine the destination rank for each migrating particle.
3. Perform the actual migration using PETSc's DMSwarmMigrate. It also calculates and logs the global number of particles migrated.

Parameters

	user	Pointer to the UserCtx structure.
	bboxlist	Array of BoundingBox structures defining the spatial domain of each MPI rank.
	migrationList_p	Pointer to a pointer for the MigrationInfo array. This array will be allocated/reallocated by IdentifyMigratingParticles if necessary. The caller is responsible for freeing this list eventually.
	migrationCount_p	Pointer to store the number of particles identified for migration on the local rank. This is reset to 0 after migration for the current cycle.
	migrationListCapacity_p	Pointer to store the current capacity of the migrationList_p array.
	currentTime	Current simulation time (used for logging).
	step	Current simulation step number (used for logging).
	migrationCycleName	A descriptive name for this migration cycle (e.g., "Preliminary Sort", "Main Loop") for logging purposes.
[out]	globalMigrationCount_out	Pointer to store the total number of particles migrated across all MPI ranks during this cycle.

Returns

PetscErrorCode 0 on success, non-zero on failure.

ReinitializeParticlesOnInletSurface()

PetscErrorCode ReinitializeParticlesOnInletSurface	(	UserCtx *	user,
		PetscReal	currentTime,
		PetscInt	step
	)

Re-initializes the positions of particles currently on this rank if this rank owns part of the designated inlet surface.

This function is intended for user->ParticleInitialization == 0 (Surface Initialization mode) and is typically called after an initial migration step (e.g., in PerformInitialSetup). It ensures that all particles that should originate from the inlet surface and are now on the correct MPI rank are properly distributed across that rank's portion of the inlet.

Parameters

user	Pointer to the UserCtx structure, containing simulation settings and grid information.
currentTime	Current simulation time (used for logging).
step	Current simulation step number (used for logging).

Returns

PetscErrorCode 0 on success, non-zero on failure.

Re-initializes the positions of particles currently on this rank if this rank owns part of the designated inlet surface.

Local to this translation unit.

Definition at line 776 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;                        // MPI rank of the current process
    DM             swarm = user->swarm;         // The particle swarm DM
    PetscReal      *positions_field = NULL;     // Pointer to swarm field for physical positions
    PetscInt64     *particleIDs = NULL;         // Pointer to swarm field for Particle IDs (for logging)
    PetscInt       *cell_ID_field = NULL;       // Pointer to swarm field for Cell IDs (for resetting after migration)
    const Cmpnts   ***coor_nodes_local_array;   // Read-only access to local node coordinates
    Vec            Coor_local;                  // Local vector for node coordinates
    DMDALocalInfo  info;                        // Local grid information (node-based) from user->da
    PetscInt       xs_gnode_rank, ys_gnode_rank, zs_gnode_rank; // Local starting node indices (incl. ghosts) of rank's DA
    PetscInt       IM_nodes_global, JM_nodes_global, KM_nodes_global; // Global node counts
 
    PetscRandom    rand_logic_reinit_i, rand_logic_reinit_j, rand_logic_reinit_k; // RNGs for re-placement
    PetscInt       nlocal_current;                // Number of particles currently on this rank
    PetscInt       particles_actually_reinitialized_count = 0; // Counter for logging
    PetscBool      can_this_rank_service_inlet = PETSC_FALSE;  // Flag
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // This function is only relevant for surface initialization mode and if an inlet face is defined.
    if ((user->simCtx->ParticleInitialization != 0 && user->simCtx->ParticleInitialization !=3) || !user->inletFaceDefined) {
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = DMSwarmGetLocalSize(swarm, &nlocal_current); CHKERRQ(ierr);
 
    // If no particles on this rank, nothing to do.
    if (nlocal_current == 0) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[T=%.4f, Step=%d] Rank %d has no local particles to re-initialize on inlet.\n", currentTime, step, rank);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    // Get DMDA information for the node-centered coordinate grid (user->da)
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    ierr = DMDAGetInfo(user->da, NULL, &IM_nodes_global, &JM_nodes_global, &KM_nodes_global, NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL); CHKERRQ(ierr);
    ierr = DMDAGetCorners(user->da, &xs_gnode_rank, &ys_gnode_rank, &zs_gnode_rank, NULL, NULL, NULL); CHKERRQ(ierr);
 
    // Modification to IM_nodes_global etc. to account for 1-cell halo in each direction.
    IM_nodes_global -= 1; JM_nodes_global -= 1; KM_nodes_global -= 1; 
 
    const PetscInt IM_cells_global = IM_nodes_global > 0 ? IM_nodes_global - 1 : 0;
    const PetscInt JM_cells_global = JM_nodes_global > 0 ? JM_nodes_global - 1 : 0;
    const PetscInt KM_cells_global = KM_nodes_global > 0 ? KM_nodes_global - 1 : 0;
 
 
 
    // Check if this rank is responsible for (part of) the designated inlet surface
    ierr = CanRankServiceInletFace(user, &info, IM_nodes_global, JM_nodes_global, KM_nodes_global, &can_this_rank_service_inlet); CHKERRQ(ierr);
 
    // Get coordinate array and swarm fields for modification
    ierr = DMGetCoordinatesLocal(user->da, &Coor_local); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, Coor_local, (void*)&coor_nodes_local_array); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),       NULL, NULL, (void**)&positions_field); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),    NULL, NULL, (void**)&particleIDs);    CHKERRQ(ierr); // For logging
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cell_ID_field);  CHKERRQ(ierr);
 
    if (!can_this_rank_service_inlet) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[T=%.4f, Step=%d] Rank %d cannot service inlet face %s. Skipping re-initialization of %d particles.\n", currentTime, step, rank, BCFaceToString(user->identifiedInletBCFace), nlocal_current);
        
        //  FALLBACK ACTION: Reset position fields to Inlet center for migration and cell ID to -1 for safety.
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[T=%.4f, Step=%d] Rank %d is resetting %d local particles to inlet center (%.6f, %.6f, %.6f) for migration.\n", currentTime, step, rank, nlocal_current, user->simCtx->CMx_c, user->simCtx->CMy_c, user->simCtx->CMz_c);
 
        for(PetscInt p = 0; p < nlocal_current; p++){
            positions_field[3*p+0] = user->simCtx->CMx_c;
            positions_field[3*p+1] = user->simCtx->CMy_c;
            positions_field[3*p+2] = user->simCtx->CMz_c;
 
            cell_ID_field[3*p+0] = -1;
            cell_ID_field[3*p+1] = -1;
            cell_ID_field[3*p+2] = -1;
        }
        
        // Cleanup: restore swarm fields/coordinate array
        ierr = DMDAVecRestoreArrayRead(user->fda, Coor_local, (void*)&coor_nodes_local_array); CHKERRQ(ierr);
        ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),       NULL, NULL, (void**)&positions_field); CHKERRQ(ierr);
        ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),    NULL, NULL, (void**)&particleIDs);    CHKERRQ(ierr); // For logging
        ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cell_ID_field);  CHKERRQ(ierr);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "[T=%.4f, Step=%d] Rank %d is on inlet face %s. Attempting to re-place %d local particles.\n", currentTime, step, rank, BCFaceToString(user->identifiedInletBCFace), nlocal_current);
 
    // Initialize fresh RNGs for this re-placement to ensure good distribution
    ierr = InitializeLogicalSpaceRNGs(user->simCtx->particleRandomSeed, &rand_logic_reinit_i, &rand_logic_reinit_j, &rand_logic_reinit_k); CHKERRQ(ierr);
 
    // Loop over all particles currently local to this rank
    for (PetscInt p = 0; p < nlocal_current; p++) {
        PetscInt  ci_metric_lnode, cj_metric_lnode, ck_metric_lnode; // Local node indices (of rank's DA patch) for cell origin
        PetscReal xi_metric_logic, eta_metric_logic, zta_metric_logic; // Intra-cell logical coordinates
        Cmpnts    phys_coords = {0.0,0.0,0.0}; // To store newly calculated physical coordinates
        PetscBool particle_was_placed = PETSC_FALSE;
 
        if(user->simCtx->ParticleInitialization == PARTICLE_INIT_SURFACE_RANDOM){
        // Get random cell on this rank's portion of the inlet and random logical coords within it
            ierr = GetRandomCellAndLogicalCoordsOnInletFace(user, &info, xs_gnode_rank, ys_gnode_rank, zs_gnode_rank,
                                                    IM_nodes_global, JM_nodes_global, KM_nodes_global,
                                                    &rand_logic_reinit_i, &rand_logic_reinit_j, &rand_logic_reinit_k,
                                                    &ci_metric_lnode, &cj_metric_lnode, &ck_metric_lnode,
                                                    &xi_metric_logic, &eta_metric_logic, &zta_metric_logic); CHKERRQ(ierr);
            
                    // Convert these logical coordinates to physical coordinates
 
            ierr = MetricLogicalToPhysical(user, coor_nodes_local_array,
                                        ci_metric_lnode, cj_metric_lnode, ck_metric_lnode,
                                        xi_metric_logic, eta_metric_logic, zta_metric_logic,
                                        &phys_coords); CHKERRQ(ierr);
                                        
            // Update the particle's position in the swarm fields
            positions_field[3*p+0] = phys_coords.x; 
            positions_field[3*p+1] = phys_coords.y; 
            positions_field[3*p+2] = phys_coords.z;
            particle_was_placed = PETSC_TRUE;                            
 
        }else if(user->simCtx->ParticleInitialization == PARTICLE_INIT_SURFACE_EDGES){
            PetscBool placement_flag = PETSC_FALSE;
            ierr = GetDeterministicFaceGridLocation(user, &info, xs_gnode_rank, ys_gnode_rank, zs_gnode_rank,
                                                IM_cells_global, JM_cells_global, KM_cells_global,
                                                particleIDs[p],
                                                &ci_metric_lnode, &cj_metric_lnode, &ck_metric_lnode,
                                                &xi_metric_logic, &eta_metric_logic, &zta_metric_logic,&placement_flag); CHKERRQ(ierr);
        
 
            if(placement_flag){
                // Convert these logical coordinates to physical coordinates
                ierr = MetricLogicalToPhysical(user, coor_nodes_local_array,
                                        ci_metric_lnode, cj_metric_lnode, ck_metric_lnode,
                                        xi_metric_logic, eta_metric_logic, zta_metric_logic,
                                        &phys_coords); CHKERRQ(ierr);
 
                // Update the particle's position in the swarm fields
                positions_field[3*p+0] = phys_coords.x;
                positions_field[3*p+1] = phys_coords.y;
                positions_field[3*p+2] = phys_coords.z;
                particle_was_placed =  PETSC_TRUE;
            } else{
                // Deterministic placement failed (particle migrated to rank where formula says it doesn't belong)
                // Fall back to random placement on this rank's portion of inlet surface
                LOG_ALLOW(GLOBAL, LOG_WARNING, "Rank %d: Particle PID %ld deterministic placement failed (belongs to different rank). Falling back to random placement.\n", rank, particleIDs[p]);
 
                ierr = GetRandomCellAndLogicalCoordsOnInletFace(user, &info, xs_gnode_rank, ys_gnode_rank, zs_gnode_rank,
                                                        IM_nodes_global, JM_nodes_global, KM_nodes_global,
                                                        &rand_logic_reinit_i, &rand_logic_reinit_j, &rand_logic_reinit_k,
                                                        &ci_metric_lnode, &cj_metric_lnode, &ck_metric_lnode,
                                                        &xi_metric_logic, &eta_metric_logic, &zta_metric_logic); CHKERRQ(ierr);
 
                // Convert to physical coordinates
                ierr = MetricLogicalToPhysical(user, coor_nodes_local_array,
                                            ci_metric_lnode, cj_metric_lnode, ck_metric_lnode,
                                            xi_metric_logic, eta_metric_logic, zta_metric_logic,
                                            &phys_coords); CHKERRQ(ierr);
 
                // Update particle position
                positions_field[3*p+0] = phys_coords.x;
                positions_field[3*p+1] = phys_coords.y;
                positions_field[3*p+2] = phys_coords.z;
                particle_was_placed = PETSC_TRUE;
            }                                    
 
        } else{
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ReinitializeParticlesOnInletSurface only supports ParticleInitialization modes 0 and 3.");
        }
 
        if(particle_was_placed){
        particles_actually_reinitialized_count++;
 
        cell_ID_field[3*p+0] = -1;
        cell_ID_field[3*p+1] = -1;
        cell_ID_field[3*p+2] = -1;
 
        LOG_LOOP_ALLOW(LOCAL, LOG_VERBOSE, p, (nlocal_current > 20 ? nlocal_current/10 : 1), // Sampled logging
            "Rank %d: PID %ld (idx %ld) RE-PLACED. CellOriginNode(locDAIdx):(%d,%d,%d). LogicCoords: (%.2e,%.2f,%.2f). PhysCoords: (%.6f,%.6f,%.6f).\n",
            rank, particleIDs[p], (long)p, 
            ci_metric_lnode, cj_metric_lnode, ck_metric_lnode,
            xi_metric_logic, eta_metric_logic, zta_metric_logic, 
            phys_coords.x, phys_coords.y, phys_coords.z);
        }
    }
 
    // Logging summary of re-initialization
    if (particles_actually_reinitialized_count > 0) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "[T=%.4f, Step=%d] Rank %d (on inlet face %d) successfully re-initialized %d of %d local particles.\n", currentTime, step, rank, user->identifiedInletBCFace, particles_actually_reinitialized_count, nlocal_current);
    } else if (nlocal_current > 0) { // This case should ideally not be hit if can_this_rank_service_inlet was true and particles were present.
         LOG_ALLOW(GLOBAL, LOG_WARNING, "[T=%.4f, Step=%d] Rank %d claimed to service inlet face %d, but re-initialized 0 of %d local particles. This may indicate an issue if particles were expected to be re-placed.\n", currentTime, step, rank, user->identifiedInletBCFace, nlocal_current);
    }
 
    // Cleanup: Destroy RNGs and restore swarm fields/coordinate array
    ierr = PetscRandomDestroy(&rand_logic_reinit_i); CHKERRQ(ierr);
    ierr = PetscRandomDestroy(&rand_logic_reinit_j); CHKERRQ(ierr);
    ierr = PetscRandomDestroy(&rand_logic_reinit_k); CHKERRQ(ierr);
 
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),       NULL, NULL, (void**)&positions_field); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),    NULL, NULL, (void**)&particleIDs);    CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID), NULL, NULL, (void**)&cell_ID_field);  CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, Coor_local, (void*)&coor_nodes_local_array); CHKERRQ(ierr);
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetLocalPIDSnapshot()

PetscErrorCode GetLocalPIDSnapshot	(	const PetscInt64	pid_field[],
		PetscInt	n_local,
		PetscInt64 **	pids_snapshot_out
	)

Creates a sorted snapshot of all Particle IDs (PIDs) from a raw data array.

This function is a crucial helper for the migration process. It captures the state of which particles are on the current MPI rank before migration occurs by taking a pointer to the swarm's raw PID data array. The resulting sorted array can then be used with an efficient binary search to quickly identify newcomer particles after migration.

This function does NOT call DMSwarmGetField/RestoreField. It is the caller's responsibility to acquire the pid_field pointer before calling and restore it afterward.

Parameters

[in]	pid_field	A read-only pointer to the raw array of PIDs for the local swarm.
[in]	n_local	The number of particles currently on the local rank.
[out]	pids_snapshot_out	A pointer to a PetscInt64* array. This function will allocate memory for this array, and the caller is responsible for freeing it with PetscFree() when it is no longer needed.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Creates a sorted snapshot of all Particle IDs (PIDs) from a raw data array.

Local to this translation unit.

Definition at line 990 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation ---
    if (!pids_snapshot_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer pids_snapshot_out is NULL.");
    }
    // If n_local > 0, pid_field must not be NULL.
    if (n_local > 0 && !pid_field) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input pid_field pointer is NULL for n_local > 0.");
    }
    
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Creating PID snapshot for %d local particles.\n", rank, n_local);
 
    // If there are no local particles, the snapshot is empty (NULL).
    if (n_local == 0) {
        *pids_snapshot_out = NULL;
 
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    
    // --- 2. Allocate Memory for the Snapshot ---
    ierr = PetscMalloc1(n_local, pids_snapshot_out); CHKERRQ(ierr);
 
    // --- 3. Copy Data ---
    // Perform a fast memory copy from the provided array to our new snapshot array.
    ierr = PetscMemcpy(*pids_snapshot_out, pid_field, n_local * sizeof(PetscInt64)); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Copied %d PIDs.\n", rank, n_local);
 
    // --- 4. Sort the Snapshot Array ---
    // Sorting enables fast binary search lookups later.
    ierr = PetscSortInt64(n_local, *pids_snapshot_out); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: PID snapshot sorted successfully.\n", rank);
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0); 
}

Here is the caller graph for this function:

AddToMigrationList()

PetscErrorCode AddToMigrationList	(	MigrationInfo **	migration_list_p,
		PetscInt *	capacity_p,
		PetscInt *	count_p,
		PetscInt	particle_local_idx,
		PetscMPIInt	destination_rank
	)

Safely adds a new migration task to a dynamically sized list.

This utility function manages a dynamic array of MigrationInfo structs. It appends a new entry to the list and automatically doubles the array's capacity using PetscRealloc if the current capacity is exceeded. This prevents buffer overflows and avoids the need to know the number of migrating particles in advance.

Parameters

[in,out]	migration_list_p	A pointer to the MigrationInfo array pointer. The function will update this pointer if the array is reallocated.
[in,out]	capacity_p	A pointer to an integer holding the current allocated capacity of the list (in number of elements). This will be updated upon reallocation.
[in,out]	count_p	A pointer to an integer holding the current number of items in the list. This will be incremented by one.
[in]	particle_local_idx	The local index (from 0 to nlocal-1) of the particle that needs to be migrated.
[in]	destination_rank	The target MPI rank for the particle.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure (e.g., from memory allocation).

Safely adds a new migration task to a dynamically sized list.

Local to this translation unit.

Definition at line 1045 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation ---
    if (!migration_list_p || !capacity_p || !count_p) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer provided to AddToMigrationList for list management.");
    }
 
    // --- 2. Check if the list needs to be resized ---
    if (*count_p >= *capacity_p) {
        PetscInt old_capacity = *capacity_p;
        // Start with a reasonable base capacity, then double for subsequent reallocations.
        PetscInt new_capacity = (old_capacity == 0) ? 16 : old_capacity * 2;
 
        // Use PetscRealloc for safe memory reallocation.
        // It handles allocating new memory, copying old data, and freeing the old block.
        // The first argument to PetscRealloc is the new size in BYTES.
        ierr = PetscRealloc(new_capacity * sizeof(MigrationInfo), migration_list_p); CHKERRQ(ierr);
        
        *capacity_p = new_capacity; // Update the capacity tracker
 
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Reallocated migrationList capacity from %d to %d.\n",
                  rank, old_capacity, new_capacity);
    }
 
    // --- 3. Add the new migration data to the list ---
    // Dereference the pointer-to-a-pointer to get the actual array.
    MigrationInfo *list = *migration_list_p;
 
    list[*count_p].local_index = particle_local_idx;
    list[*count_p].target_rank = destination_rank;
 
    // --- 4. Increment the count of items in the list ---
    (*count_p)++;
    
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

FlagNewcomersForLocation()

PetscErrorCode FlagNewcomersForLocation	(	DM	swarm,
		PetscInt	n_local_before,
		const PetscInt64	pids_before[]
	)

Identifies newly arrived particles after migration and flags them for a location search.

This function is a critical component of the iterative migration process managed by the main particle settlement orchestrator (e.g., SettleParticles). After a DMSwarmMigrate call, each rank's local particle list is a new mix of resident particles and newly received ones. This function's job is to efficiently identify these "newcomers" and set their DMSwarm_location_status field to NEEDS_LOCATION.

This ensures that in the subsequent pass of the migration do-while loop, only the newly arrived particles are processed by the expensive location algorithm, preventing redundant work on particles that are already settled on the current rank.

The identification is done by comparing the PIDs of particles currently on the rank against a "snapshot" of PIDs taken before the migration occurred.

Parameters

[in]	swarm	The DMSwarm object, which has just completed a migration.
[in]	n_local_before	The number of particles that were on this rank before the migration was performed.
[in]	pids_before	A pre-sorted array of the PIDs that were on this rank before the migration. This is used for fast lookups.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Note

This function assumes the pids_before array is sorted in ascending order to enable the use of an efficient binary search.

Identifies newly arrived particles after migration and flags them for a location search.

Local to this translation unit.

Definition at line 1103 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscInt       n_local_after;
    PetscInt       newcomer_count = 0;
    
    // Pointers to the swarm data fields we will read and modify
    PetscInt64 *pid_field_after    = NULL;
    PetscInt   *status_field_after = NULL;
    PetscInt   *cell_field_after = NULL;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation and Basic Setup ---
    if (!swarm) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input DMSwarm is NULL in FlagNewcomersForLocation.");
    }
    // If n_local_before > 0, the corresponding PID array must not be null.
    if (n_local_before > 0 && !pids_before) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input pids_before array is NULL for n_local_before > 0.");
    }
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    
    // Get the number of particles on this rank *after* the migration.
    ierr = DMSwarmGetLocalSize(swarm, &n_local_after); CHKERRQ(ierr);
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Checking for newcomers. Size before: %d, Size after: %d\n",
              rank, n_local_before, n_local_after);
 
    // If there are no particles now, there's nothing to do.
    if (n_local_after == 0) {
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    
    // --- 2. Access Swarm Data ---
    // Get read-only access to the PIDs and read-write access to the status field.
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_field_after); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_field_after); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_field_after); CHKERRQ(ierr);
    if (!pid_field_after || !status_field_after) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Failed to get required swarm fields in FlagNewcomersForLocation.");
    }
 
    // --- 3. Identify and Flag Newcomers ---
    // Loop through all particles currently on this rank.
    for (PetscInt p_idx = 0; p_idx < n_local_after; ++p_idx) {
        PetscInt64 current_pid = pid_field_after[p_idx];
        PetscBool  is_found_in_before_list;
 
        // Use our custom, efficient helper function for the lookup.
        ierr = BinarySearchInt64(n_local_before, pids_before, current_pid, &is_found_in_before_list); CHKERRQ(ierr);
 
        // If the PID was NOT found in the "before" list, it must be a newcomer.
        if (!is_found_in_before_list) {
          // Flag it for processing in the next pass of the migration loop.
          status_field_after[p_idx] = NEEDS_LOCATION;
      // cell_field_after[3*p_idx+0] = -1;
      // cell_field_after[3*p_idx+1] = -1; 
      // cell_field_after[3*p_idx+2] = -1; 
            newcomer_count++;
            
            LOG_ALLOW(LOCAL, LOG_VERBOSE, "[Rank %d]: Flagged newcomer PID %ld at local index %d as NEEDS_LOCATION.\n",
                      rank, current_pid, p_idx);
        }
    }
 
    // --- 4. Restore Swarm Fields ---
    // Release the locks on the swarm data arrays.
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_field_after); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_field_after); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_field_after); CHKERRQ(ierr);
 
    if (newcomer_count > 0) {
        LOG_ALLOW(LOCAL, LOG_INFO, "[Rank %d]: Identified and flagged %d newcomers.\n", rank, newcomer_count);
    }
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MigrateRestartParticlesUsingCellID()

PetscErrorCode MigrateRestartParticlesUsingCellID

(

UserCtx *

user

)

Fast-path migration for restart particles using preloaded Cell IDs.

This function provides an optimized migration path specifically for particles loaded from restart files. Unlike the standard LocateAllParticlesInGrid() which performs expensive walking searches, this function leverages the fact that restart particles already have valid global Cell IDs loaded from disk.

How It Works:

1. Iterates through all local particles.
2. For each particle with a valid Cell ID (ci, cj, ck):
   • Calls FindOwnerOfCell(ci, cj, ck) to determine the correct rank.
   • If owner differs from current rank, adds to migration list.
   • If owner matches current rank, the existing ACTIVE_AND_LOCATED status is preserved.
3. Uses existing SetMigrationRanks() and PerformMigration() infrastructure.
4. Achieves single-pass direct migration (no multi-hop, no walking searches).

Parameters

[in,out]

user

Pointer to UserCtx containing the swarm and RankCellInfoMap. The function updates particle status fields and performs migration.

Returns

PetscErrorCode 0 on success, non-zero on failure.

Note

Testing status: Direct coverage currently focuses on restart fast-path ownership transfer. Non-restart multi-pass migration behavior remains part of the next simulation-core test backlog.

Fast-path migration for restart particles using preloaded Cell IDs.

Local to this translation unit.

Definition at line 1197 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    PetscInt       nlocal;
    PetscInt       *cell_p = NULL;
    PetscInt64     *pid_p = NULL;
    PetscMPIInt    rank;
    
    MigrationInfo  *migrationList = NULL;
    PetscInt       local_migration_count = 0;
    PetscInt       migrationListCapacity = 0;
    PetscInt       global_migration_count = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Checking %d restart particles for direct migration using CellIDs.\n", nlocal);
 
    if (nlocal > 0) {
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_p);   CHKERRQ(ierr);
        ierr = DMSwarmGetField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_p);    CHKERRQ(ierr);
        
        // Note: We do NOT need to modify the status field here. 
        // We trust the loaded status (ACTIVE_AND_LOCATED) is correct for the destination rank.
 
        for (PetscInt p_idx = 0; p_idx < nlocal; ++p_idx) {
            PetscInt ci = cell_p[3*p_idx + 0];
            PetscInt cj = cell_p[3*p_idx + 1];
            PetscInt ck = cell_p[3*p_idx + 2];
 
            /* Skip particles with invalid Cell IDs (will be handled by LocateAllParticles) */
            if (ci < 0 || cj < 0 || ck < 0) {
                continue;
            }
 
            PetscMPIInt owner_rank;
            ierr = FindOwnerOfCell(user, ci, cj, ck, &owner_rank); CHKERRQ(ierr);
 
            if (owner_rank != -1 && owner_rank != rank) {
                /* Particle belongs to another rank - migrate it */
                ierr = AddToMigrationList(&migrationList, &migrationListCapacity, &local_migration_count,
                                          p_idx, owner_rank); CHKERRQ(ierr);
                
                LOG_ALLOW(LOCAL, LOG_VERBOSE, "[PID %ld] Direct migration: Cell (%d,%d,%d) belongs to Rank %d (Current: %d).\n",
                          (long)pid_p[p_idx], ci, cj, ck, owner_rank, rank);
            }
        }
 
        ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_p);   CHKERRQ(ierr);
        ierr = DMSwarmRestoreField(swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_p);    CHKERRQ(ierr);
    }
 
    /* Check if any rank needs to migrate particles */
    ierr = MPI_Allreduce(&local_migration_count, &global_migration_count, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    if (global_migration_count > 0) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Fast restart migration: Directly migrating %d particles using CellIDs.\n", global_migration_count);
        ierr = SetMigrationRanks(user, migrationList, local_migration_count); CHKERRQ(ierr);
        ierr = PerformMigration(user); CHKERRQ(ierr);
        /* We do NOT flag newcomers here. We trust their loaded status (ACTIVE_AND_LOCATED) */
        /* is valid for their destination rank. */
    } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Fast restart migration: All particles are already on correct ranks.\n");
    }
 
    ierr = PetscFree(migrationList); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LocateAllParticlesInGrid()

PetscErrorCode LocateAllParticlesInGrid	(	UserCtx *	user,
		BoundingBox *	bboxlist
	)

Orchestrates the complete particle location and migration process for one timestep.

This function is the master orchestrator for ensuring every particle is on its correct MPI rank and has a valid host cell index. It is designed to be called once per timestep after particle positions have been updated.

The function uses a robust, iterative "Guess and Verify" strategy within a do-while loop to handle complex particle motion across processor boundaries, especially on curvilinear grids.

1. State Snapshot: At the start of each pass, it captures a list of all Particle IDs (PIDs) on the current rank.
2. **"Guess" (Heuristic):** For particles that are "lost" (no valid host cell), it first attempts a fast, bounding-box-based guess to find a potential new owner rank.
3. **"Verify" (Robust Walk):** For all other particles, or if the guess fails, it uses a robust cell-walking algorithm (LocateParticleOrFindMigrationTarget) that determines the particle's status: located locally, needs migration, or is lost.
4. Migration: After identifying all migrating particles on a pass, it performs the MPI communication using the SetMigrationRanks and PerformMigration helpers.
5. Newcomer Flagging: After migration, it uses the PID snapshot from step 1 to efficiently identify newly arrived particles and flag them for location on the next pass.
6. Iteration: The process repeats in a do-while loop until a pass occurs where no particles migrate, ensuring the entire swarm is in a stable, consistent state.

Parameters

[in,out]	user	Pointer to the UserCtx, containing the swarm and all necessary domain topology information (bboxlist, RankCellInfoMap, etc.).
[in]	bboxlist	An array of BoundingBox structures for ALL MPI ranks, indexed 0 to (size-1). This array must be up-to-date and available on all ranks.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Note

Testing status: Direct unit coverage currently pins the prior-cell fast path and the local guess-then-verify path. Multi-pass migration, newcomer flagging, and several lost/migration edge cases are still targeted for future bespoke tests.

Orchestrates the complete particle location and migration process for one timestep.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/ParticleMotion.h.

See also

LocateAllParticlesInGrid()

Definition at line 1399 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscInt       passes = 0;
    const PetscInt MAX_MIGRATION_PASSES = 50; // Safety break for runaway loops
    PetscInt       global_migrations_this_pass;
    PetscMPIInt    rank;
    PetscInt total_migrated_this_timestep = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = ResetSearchMetrics(user->simCtx); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "LocateAllParticlesInGrid (Orchestrator) - Beginning particle settlement process.\n");
 
    // This loop ensures that particles that jump across multiple ranks are
    // handled correctly in successive, iterative handoffs.
    do {
        passes++;
        user->simCtx->searchMetrics.currentSettlementPass = passes;
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Starting migration pass %d.\n", rank, passes);
 
        // --- STAGE 1: PER-PASS INITIALIZATION ---
        MigrationInfo  *migrationList = NULL;
        PetscInt       local_migration_count = 0;
        PetscInt       migrationListCapacity = 0;
        PetscInt       nlocal_before;
        PetscInt64     *pids_before_snapshot = NULL;
        PetscInt       local_lost_count = 0;
 
        ierr = DMSwarmGetLocalSize(user->swarm, &nlocal_before); CHKERRQ(ierr);
        if (passes == 1) {
            user->simCtx->searchMetrics.searchPopulation += (PetscInt64)nlocal_before;
        }
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Pass %d begins with %d local particles.\n", rank, passes, nlocal_before);
 
 
        // --- STAGE 2: PRE-MIGRATION SNAPSHOT & MAIN PROCESSING LOOP ---
        if (nlocal_before > 0) {
            // Get pointers to all fields needed for this pass
            PetscReal  *pos_p, *weights_p, *vel_p;
            PetscInt   *cell_p, *status_p;
            PetscInt64 *pid_p;
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void**)&pos_p);    CHKERRQ(ierr);
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY),                NULL, NULL, (void**)&vel_p);    CHKERRQ(ierr);
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT),                  NULL, NULL, (void**)&weights_p);  CHKERRQ(ierr);
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_p);     CHKERRQ(ierr);
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_p);      CHKERRQ(ierr);
            ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_p);   CHKERRQ(ierr);
 
            // Create a sorted snapshot of current PIDs to identify newcomers after migration.
            // This helper requires a raw pointer, which we just acquired.
            ierr = GetLocalPIDSnapshot(pid_p, nlocal_before, &pids_before_snapshot); CHKERRQ(ierr);
 
            for (PetscInt p_idx = 0; p_idx < nlocal_before; p_idx++) {
 
                // OPTIMIZATION: Skip particles already settled in a previous pass of this do-while loop.
 
          LOG_ALLOW(LOCAL,LOG_VERBOSE,
            "Local Particle idx=%d, PID=%ld, status=%s, cell=(%d, %d, %d)\n",
            p_idx,
            (long)pid_p[p_idx],
            ParticleLocationStatusToString((ParticleLocationStatus)status_p[p_idx]),
            cell_p[3*p_idx],
            cell_p[3*p_idx+1],
            cell_p[3*p_idx+2]);
 
                if (status_p[p_idx] == ACTIVE_AND_LOCATED) {
          LOG_ALLOW(LOCAL,LOG_VERBOSE," [rank %d][PID %ld] skipped in pass %d as it is already located at (%d,%d,%d).\n",rank,pid_p[p_idx],passes,cell_p[3*p_idx],cell_p[3*p_idx + 1],cell_p[3*p_idx + 2]);
                    continue;
                }
 
                // UNPACK: Create a temporary C struct for easier processing using our helper.
                Particle current_particle;
 
            //  LOG_ALLOW(LOCAL,LOG_DEBUG,"about to unpack p_idx=%d (PID=%ld)\n",p_idx, (long)pid_p[p_idx]);
        
                ierr = UnpackSwarmFields(p_idx, pid_p, weights_p, pos_p, cell_p, vel_p, status_p,NULL,NULL,NULL,&current_particle); CHKERRQ(ierr);
 
            //      LOG_ALLOW(LOCAL,LOG_DEBUG,"unpacked p_idx=%d → cell[0]=%d, status=%s\n",p_idx, current_particle.cell[0], ParticleLocationStatusToString((ParticleLocationStatus)current_particle.location_status));
 
            ParticleLocationStatus final_status = (ParticleLocationStatus)status_p[p_idx];
 
 
            // CASE 1: Particle has a valid prior cell index.
                // It has moved, so we only need to run the robust walk from its last known location.
                if (current_particle.cell[0] >= 0) {
                    LOG_ALLOW(LOCAL, LOG_VERBOSE, "[PID %ld] has valid prior cell. Strategy: Robust Walk from previous cell.\n", current_particle.PID);
                    ierr = LocateParticleOrFindMigrationTarget(user, &current_particle, &final_status); CHKERRQ(ierr);
                } 
 
        /*      
                // --- "GUESS" FAST PATH for lost particles ---
                if (current_particle.cell[0] < 0) {
                    LOG_ALLOW(LOCAL, LOG_DEBUG, "[PID %ld] is lost or uninitialzied (cell=%d), attempting fast guess.\n",current_particle.PID, current_particle.cell[0]);
                    ierr = GuessParticleOwnerWithBBox(user, &current_particle, bboxlist, &destination_rank); CHKERRQ(ierr);
                    if (destination_rank != MPI_PROC_NULL && destination_rank != rank) {
                        final_status = MIGRATING_OUT;
                        // The particle struct's destination rank must be updated for consistency
                        current_particle.destination_rank = destination_rank;
                    }
                }
 
        LOG_ALLOW(LOCAL,LOG_DEBUG,"[PID %ld] Particle status after Initial Guess:%d \n",current_particle.PID,final_status);
 
                // --- "VERIFY" ROBUST WALK if guess didn't resolve it ---
                if (final_status == NEEDS_LOCATION  || UNINITIALIZED) {
                    LOG_ALLOW(LOCAL, LOG_DEBUG, "[PID %ld] Not resolved by guess, starting robust walk.\n", current_particle.PID);
                    // This function will update the particle's status and destination rank internally.
             ierr = LocateParticleOrFindMigrationTarget(user, &current_particle, &final_status); CHKERRQ(ierr);
                    destination_rank = current_particle.destination_rank; // Retrieve the result
                }
 
                // --- PROCESS THE FINAL STATUS AND TAKE ACTION ---
                if (final_status == MIGRATING_OUT) {
                    status_p[p_idx] = MIGRATING_OUT; // Mark for removal by DMSwarm
                    ierr = AddToMigrationList(&migrationList, &migrationListCapacity, &local_migration_count, p_idx, destination_rank); CHKERRQ(ierr);
                    LOG_ALLOW(LOCAL, LOG_DEBUG, "[PID %ld] at local index %d marked for migration to rank %d.\n",current_particle.PID, p_idx, destination_rank);
                } else {
                     // Particle's final status is either LOCATED or LOST; update its state in the swarm arrays.
                     current_particle.location_status = final_status;
                     // PACK: Use the helper to write results back to the swarm arrays.
                     ierr = UpdateSwarmFields(p_idx, &current_particle, pos_p, vel_p, weights_p, cell_p, status_p,NULL,NULL,NULL); CHKERRQ(ierr);
                }
        */
                // CASE 2: Particle is "lost" (cell = -1). Strategy: Guess -> Verify.
                else {
          LOG_ALLOW(LOCAL, LOG_VERBOSE, "[PID %ld] has invalid cell. Strategy: Guess Owner -> Find Cell.\n",current_particle.PID);
                    
          PetscMPIInt guessed_owner_rank = MPI_PROC_NULL;
          ierr = GuessParticleOwnerWithBBox(user, &current_particle, bboxlist, &guessed_owner_rank); CHKERRQ(ierr);
 
          // If the guess finds a DIFFERENT rank, we can mark for migration and skip the walk.
          if (guessed_owner_rank != MPI_PROC_NULL && guessed_owner_rank != rank) {
            user->simCtx->searchMetrics.bboxGuessSuccessCount++;
            LOG_ALLOW(LOCAL, LOG_VERBOSE, "[PID %ld] Guess SUCCESS: Found migration target Rank %d. Finalizing.\n", current_particle.PID, guessed_owner_rank);
            final_status = MIGRATING_OUT;
            current_particle.destination_rank = guessed_owner_rank;
          } 
          else {
            user->simCtx->searchMetrics.bboxGuessFallbackCount++;
 
            // This block runs if the guess either failed (rank is NULL) or found the particle is local (rank is self).
            // In BOTH cases, the situation is unresolved, and we MUST fall back to the robust walk.
            if (guessed_owner_rank == rank) {
              LOG_ALLOW(LOCAL, LOG_DEBUG, "[PID %ld] Guess determined particle is local. Proceeding to robust walk to find cell.\n", current_particle.PID);
            } else { // guessed_owner_rank == MPI_PROC_NULL
              LOG_ALLOW(LOCAL, LOG_WARNING, "[PID %ld] Guess FAILED to find an owner. Proceeding to robust walk for definitive search.\n", current_particle.PID);
            }
                        
            ierr = LocateParticleOrFindMigrationTarget(user, &current_particle, &final_status); CHKERRQ(ierr);
          }
                }
        
                // --- PROCESS THE FINAL, DEFINITIVE STATUS ---
                current_particle.location_status = final_status;
                ierr = UpdateSwarmFields(p_idx, &current_particle, pos_p, vel_p, weights_p, cell_p, status_p,NULL,NULL,NULL); CHKERRQ(ierr);
                
                if (final_status == MIGRATING_OUT) {
          ierr = AddToMigrationList(&migrationList, &migrationListCapacity, &local_migration_count, p_idx, current_particle.destination_rank); CHKERRQ(ierr);
                } else if (final_status == LOST) {
          local_lost_count++;
                  user->simCtx->searchMetrics.searchLostCount++;
                } else if (final_status == ACTIVE_AND_LOCATED) {
                  user->simCtx->searchMetrics.searchLocatedCount++;
                }
                
            } // End of main particle processing loop
 
            // Restore all the fields acquired for this pass.
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),                NULL, NULL, (void**)&pos_p);    CHKERRQ(ierr);
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_VELOCITY),                NULL, NULL, (void**)&vel_p);    CHKERRQ(ierr);
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_WEIGHT),                  NULL, NULL, (void**)&weights_p);  CHKERRQ(ierr);
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_CELL_ID),          NULL, NULL, (void**)&cell_p);     CHKERRQ(ierr);
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_PID),             NULL, NULL, (void**)&pid_p);      CHKERRQ(ierr);
            ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_p);   CHKERRQ(ierr);
        }
 
        // --- STAGE 3: ACTION & MPI COMMUNICATION ---
        LOG_ALLOW(LOCAL, LOG_INFO, "[Rank %d] Pass %d: Identified %d particles to migrate out.\n", rank, passes, local_migration_count);
 
        // --- STAGE 3: SYNCHRONIZE AND DECIDE ---
        // FIRST, determine if any rank wants to migrate. This call is safe because
        // all ranks have finished their local work and can participate.
        ierr = MPI_Allreduce(&local_migration_count, &global_migrations_this_pass, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);  
 
        total_migrated_this_timestep += global_migrations_this_pass;
 
        if(global_migrations_this_pass > 0 ){
 
            LOG_ALLOW(GLOBAL, LOG_INFO, "Pass %d: Migrating %d particles globally.\n", passes, global_migrations_this_pass);
                
                ierr = SetMigrationRanks(user, migrationList, local_migration_count); CHKERRQ(ierr);
                ierr = PerformMigration(user); CHKERRQ(ierr);
 
                // --- STAGE 4: POST-MIGRATION RESET ---
                // Identify newly arrived particles and flag them with NEEDS_LOCATION so they are
                // processed in the next pass. This uses the snapshot taken in STAGE 2.
                ierr = FlagNewcomersForLocation(user->swarm, nlocal_before, pids_before_snapshot); CHKERRQ(ierr);
        }
        // --- STAGE 5: LOOP SYNCHRONIZATION AND CLEANUP ---
 
        ierr = PetscFree(pids_before_snapshot);
        ierr = PetscFree(migrationList);
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "End of pass %d. Total particles migrated globally: %d.\n", passes, global_migrations_this_pass);
 
    } while (global_migrations_this_pass > 0 && passes < MAX_MIGRATION_PASSES);
 
    // --- FINAL CHECKS ---
    if (passes >= MAX_MIGRATION_PASSES) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_CONV_FAILED, "Particle migration failed to converge after %d passes. Check for particles oscillating between ranks.", MAX_MIGRATION_PASSES);
    }
 
    user->simCtx->particlesMigratedLastStep = total_migrated_this_timestep;
    user->simCtx->migrationPassesLastStep = passes;
    user->simCtx->searchMetrics.maxParticlePassDepth = PetscMax(user->simCtx->searchMetrics.maxParticlePassDepth, (PetscInt64)passes);
    user->simCtx->searchMetrics.currentSettlementPass = 0;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Location completed in %d passes.\n", passes);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ResetAllParticleStatuses()

PetscErrorCode ResetAllParticleStatuses

(

UserCtx *

user

)

Marks all local particles as NEEDS_LOCATION for the next settlement pass.

This function is designed to be called at the end of a full timestep, after all particle-based calculations are complete. It prepares the swarm for the next timestep by ensuring that after the next position update, every particle will be re-evaluated by the LocateAllParticlesInGrid orchestrator.

It iterates through all locally owned particles and sets their DMSwarm_location_status field to NEEDS_LOCATION.

Parameters

[in,out]

user

Pointer to the UserCtx containing the swarm.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Marks all local particles as NEEDS_LOCATION for the next settlement pass.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/ParticleMotion.h.

See also

ResetAllParticleStatuses()

Definition at line 1633 of file ParticleMotion.c.

{
    PetscErrorCode ierr;
    PetscInt       n_local;
    PetscInt      *status_p;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
 
    if (n_local > 0) {
        // Get write access to the status field
        ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_p); CHKERRQ(ierr);
        
        for (PetscInt p = 0; p < n_local; ++p) {
            // Only reset particles that are considered settled. This is a small optimization
            // to avoid changing the status of a LOST particle, though resetting all would also be fine.
            if (status_p[p] == ACTIVE_AND_LOCATED) {
                status_p[p] = NEEDS_LOCATION;
            }
        }
        
        // Restore the field
        ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_LOCATION_STATUS), NULL, NULL, (void**)&status_p); CHKERRQ(ierr);
    }
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: