Skip to content

Commit

Permalink
Particle destructor (#927)
Browse files Browse the repository at this point in the history
### Description
Support particle destruction. 

Changes:
- Define `ParticleDestructionTraits`. The user define their own particle
destructor by specializing the `ParticleChecker` method in
`ParticleDestructionTraits`, which return true or false to mark a
particle for destruction. As a placeholder, `particles.param3` is used
in particle destructor.
- Define `ParticleDestructionImpl`, which loops over all particles at a
level and mark particles for destruction.
- Add a `getNumParticles` function to `PhysicsParticleDescriptor` to get
the total number of *valid* particles of that type.
- Add a new particle parameter `particles.verbose` to turn on printing
particle logistics at each time step. Currently, it prints the number of
particles of each type.

The new particle destructor feature is tested in gravity_3d.cpp. 

Next, I'll add a `particleDestruction` method in
`ParticleDestructionTraits`, which modifies the hydro state while
removing a particle.

### Related issues
None. 

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues (if applicable; see
above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [ ] *(For quokka-astro org members)* I have manually triggered the GPU
tests with the magic comment `/azp run`.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
  • Loading branch information
chongchonghe and pre-commit-ci[bot] authored Mar 7, 2025
1 parent 824c6c7 commit 306cf35
Show file tree
Hide file tree
Showing 7 changed files with 219 additions and 23 deletions.
66 changes: 56 additions & 10 deletions src/particles/PhysicsParticles.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "AMReX_Vector.H"
#include "particle_creation.hpp"
#include "particle_deposition.hpp"
#include "particle_destruction.hpp"
#include "particle_types.hpp"
#include "physics_info.hpp"

Expand All @@ -32,10 +33,12 @@ class PhysicsParticleDescriptorBase
int birthTimeIndex_{-1}; // Index for birth time (-1 if not used)
bool interactsWithHydro_{false}; // Whether particles interact with hydrodynamics
bool allowsCreation_{false}; // Whether particles can be created during simulation
bool allowsDestruction_{false}; // Whether particles can be destroyed during simulation

public:
PhysicsParticleDescriptorBase(int mass_idx, int lum_idx, int birth_time_idx, bool hydro_interact, bool allows_creation)
: massIndex_(mass_idx), lumIndex_(lum_idx), birthTimeIndex_(birth_time_idx), interactsWithHydro_(hydro_interact), allowsCreation_(allows_creation)
PhysicsParticleDescriptorBase(int mass_idx, int lum_idx, int birth_time_idx, bool hydro_interact, bool allows_creation, bool allows_destruction = false)
: massIndex_(mass_idx), lumIndex_(lum_idx), birthTimeIndex_(birth_time_idx), interactsWithHydro_(hydro_interact), allowsCreation_(allows_creation),
allowsDestruction_(allows_destruction)
{
}

Expand All @@ -53,6 +56,7 @@ class PhysicsParticleDescriptorBase
[[nodiscard]] AMREX_FORCE_INLINE auto getBirthTimeIndex() const -> int { return birthTimeIndex_; }
[[nodiscard]] AMREX_FORCE_INLINE auto getInteractsWithHydro() const -> bool { return interactsWithHydro_; }
[[nodiscard]] AMREX_FORCE_INLINE auto getAllowsCreation() const -> bool { return allowsCreation_; }
[[nodiscard]] AMREX_FORCE_INLINE auto getAllowsDestruction() const -> bool { return allowsDestruction_; }

// Virtual interface for particle operations
[[nodiscard]] virtual auto getParticlePositions(int lev) const -> std::vector<std::array<double, AMREX_SPACEDIM>> = 0;
Expand All @@ -66,11 +70,13 @@ class PhysicsParticleDescriptorBase
virtual void redistribute(int lev, int ngrow) = 0;
virtual void writePlotFile(const std::string &plotfilename, const std::string &name) = 0;
virtual void writeCheckpoint(const std::string &checkpointname, const std::string &name, bool include_header) = 0;
[[nodiscard]] virtual auto getNumParticles() const -> int = 0;
#if AMREX_SPACEDIM == 3
virtual void depositMass(const amrex::Vector<amrex::MultiFab *> &rhs, int finest_lev, amrex::Real Gconst) = 0;
virtual void driftParticles(int lev, amrex::Real dt) const = 0;
virtual void kickParticles(int lev, amrex::Real dt, amrex::MultiFab const &acceleration) = 0;
virtual void createParticlesFromState(amrex::MultiFab &state, int lev, amrex::Real current_time, amrex::Real dt) const = 0;
virtual void destroyParticles(int lev, amrex::Real current_time, amrex::Real dt) = 0;
[[nodiscard]] virtual auto computeMaxParticleSpeed(int lev) const -> amrex::Real = 0;
#endif // AMREX_SPACEDIM == 3
};
Expand All @@ -87,8 +93,9 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
[[nodiscard]] static constexpr auto getParticleType() -> ParticleType { return particleType_; }

// Constructor initializing descriptor with container and particle properties
PhysicsParticleDescriptor(int mass_idx, int lum_idx, int birth_time_idx, bool hydro_interact, bool allows_creation, ContainerType *container)
: PhysicsParticleDescriptorBase(mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation), container_(container)
PhysicsParticleDescriptor(int mass_idx, int lum_idx, int birth_time_idx, bool hydro_interact, bool allows_creation, ContainerType *container,
bool allows_destruction = false)
: PhysicsParticleDescriptorBase(mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, allows_destruction), container_(container)
{
}

Expand Down Expand Up @@ -209,6 +216,15 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
return particle_data; // Empty vector on non-root ranks
}

// Get the number of particles in the container
[[nodiscard]] auto getNumParticles() const -> int override
{
if (container_ != nullptr) {
return static_cast<int>(container_->TotalNumberOfParticles(true, false));
}
return 0;
}

#if AMREX_SPACEDIM == 3

// Implementation of mass deposition from particles to grid
Expand Down Expand Up @@ -293,6 +309,14 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
current_time, dt);
}

void destroyParticles(int lev, amrex::Real current_time, amrex::Real dt) override
{
if (container_ != nullptr) {
ParticleDestructionTraits<particleType_>::template destroyParticles<problem_t, ContainerType>(container_, this->getMassIndex(), lev,
current_time, dt);
}
}

// Compute maximum particle speed at a given level
[[nodiscard]] auto computeMaxParticleSpeed(int lev) const -> amrex::Real override
{
Expand Down Expand Up @@ -422,22 +446,22 @@ template <typename problem_t> class PhysicsParticleRegister
// Register a new particle type with specified properties
template <typename ContainerType>
void registerParticleType(ParticleType type, int mass_idx, int lum_idx, int birth_time_idx, bool hydro_interact, bool allows_creation,
ContainerType *container)
ContainerType *container, bool allows_destruction = false)
{
std::unique_ptr<PhysicsParticleDescriptorBase> descriptor;

// Create the appropriate descriptor based on the particle type
if (type == ParticleType::Rad) {
descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::Rad>>(
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container);
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container, allows_destruction);
}
#if AMREX_SPACEDIM == 3
else if (type == ParticleType::CIC) {
descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::CIC>>(
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container);
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container, allows_destruction);
} else if (type == ParticleType::CICRad) {
descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::CICRad>>(
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container);
mass_idx, lum_idx, birth_time_idx, hydro_interact, allows_creation, container, allows_destruction);
}
#endif // AMREX_SPACEDIM == 3
else {
Expand Down Expand Up @@ -547,11 +571,23 @@ template <typename problem_t> class PhysicsParticleRegister
}
}

// Compute maximum particle speed at a given level
// Destroy particles based on particle type
void destroyParticles(int lev, amrex::Real current_time, amrex::Real dt)
{
for (const auto &[type, descriptor] : particleRegistry_) {
// Only destroy particles if the descriptor allows destruction
if (descriptor->getAllowsDestruction()) {
// Call the appropriate particle destruction method based on the particle type
descriptor->destroyParticles(lev, current_time, dt);
}
}
}

// Compute maximum particle speed across all particle types
[[nodiscard]] auto computeMaxParticleSpeed(int lev) const -> amrex::Real
{
amrex::Real max_speed = 0.0;
for (const auto &[name, descriptor] : particleRegistry_) {
for (const auto &[type, descriptor] : particleRegistry_) {
if (descriptor->getMassIndex() >= 0) {
const amrex::Real speed = descriptor->computeMaxParticleSpeed(lev);
AMREX_ASSERT(!std::isnan(speed));
Expand All @@ -562,6 +598,16 @@ template <typename problem_t> class PhysicsParticleRegister
}
#endif // AMREX_SPACEDIM == 3

// Print particle statistics
void printParticleStatistics() const
{
amrex::Print() << "Particle statistics:\n";
amrex::Print() << "Particle type, Number of particles\n";
for (const auto &[type, descriptor] : particleRegistry_) {
amrex::Print() << getParticleTypeName(type) << ", " << descriptor->getNumParticles() << "\n";
}
}

// Prevent copying or moving of the registry to ensure single ownership
PhysicsParticleRegister(const PhysicsParticleRegister &) = delete;
auto operator=(const PhysicsParticleRegister &) -> PhysicsParticleRegister & = delete;
Expand Down
2 changes: 1 addition & 1 deletion src/particles/particle_creation.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ template <typename problem_t, typename ContainerType, template <typename> class
static void createParticlesImpl(ContainerType *container, int mass_idx, amrex::MultiFab &state, int lev, amrex::Real current_time, amrex::Real dt)
{
if (container != nullptr) {
if (mass_idx >= 0 && mass_idx + 3 < ContainerType::ParticleType::NReal) {
if (mass_idx >= 0) {
// Use the provided ParticleChecker type with global particle parameters
CheckerType<problem_t> particle_checker;

Expand Down
82 changes: 82 additions & 0 deletions src/particles/particle_destruction.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
#ifndef PARTICLE_DESTRUCTION_HPP_
#define PARTICLE_DESTRUCTION_HPP_

#include "particle_types.hpp"

namespace quokka
{

// Helper namespace with implementation details for particle destruction
namespace ParticleDestructionImpl
{
// Common implementation of particle destruction logic
template <typename problem_t, typename ContainerType, template <typename> class CheckerType>
static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev, amrex::Real current_time, amrex::Real dt)
{
if (container != nullptr) {
if (mass_idx >= 0) {
// Use the provided ParticleChecker type to determine which particles to destroy
CheckerType<problem_t> particle_checker;

// Iterate through all particles at this level
for (typename ContainerType::ParIterType pti(*container, lev); pti.isValid(); ++pti) {
auto &particles = pti.GetArrayOfStructs();
const int np = particles.numParticles();

// Skip if no particles in this tile
if (np == 0) {
continue;
}

// Process particles on the device
auto *parray = particles().data();

amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) {
auto &p = parray[i]; // NOLINT

// Check if this particle should be destroyed
if (particle_checker(p, mass_idx, current_time, dt)) {
// Mark particle as invalid by negating its ID
// This is the AMReX way to mark particles for removal
// They will be removed during the next Redistribute call
p.id() = -1;
}
});
}

// Redistribute particles at this level to actually remove the invalid particles
// TODO(cch): This won't work when AMR subcycling is enabled. When a particle moves into the ghost cells in the first step of the
// subcycle, it may be moved from that level into a lower level. Then, in the second step of the subcycle, it will not be drifted.
container->Redistribute(lev);
}
}
}
} // namespace ParticleDestructionImpl

// Traits class for specializing particle destruction behavior
template <ParticleType particleType> struct ParticleDestructionTraits {
// Default nested ParticleChecker - determines if a particle should be destroyed
template <typename problem_t> struct ParticleChecker {

template <typename ParticleType>
AMREX_GPU_DEVICE auto operator()(ParticleType &p, int mass_idx, amrex::Real current_time, amrex::Real dt) const -> bool
{
// Default implementation: destroy particles with mass < 1.0
amrex::ignore_unused(p, mass_idx, current_time, dt);
return false;
}
};

// Main method to destroy particles - uses the helper implementation
template <typename problem_t, typename ContainerType>
static void destroyParticles(ContainerType *container, int mass_idx, int lev, amrex::Real current_time, amrex::Real dt)
{
// Use the common implementation with our checker type
ParticleDestructionImpl::destroyParticlesImpl<problem_t, ContainerType, ParticleDestructionTraits<particleType>::template ParticleChecker>(
container, mass_idx, lev, current_time, dt);
}
};

} // namespace quokka

#endif // PARTICLE_DESTRUCTION_HPP_
4 changes: 4 additions & 0 deletions src/particles/particle_types.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,8 @@ enum class ParticleType {
// their own copies.
inline amrex::Real particle_param1 = -1.0; // NOLINT
inline amrex::Real particle_param2 = -1.0; // NOLINT
inline amrex::Real particle_param3 = -1.0; // NOLINT
inline int particle_verbose = 0; // NOLINT print particle logistics

//-------------------- Radiation particles --------------------

Expand Down Expand Up @@ -166,6 +168,8 @@ inline void particleParmParse()
const amrex::ParmParse pp("particles");
pp.query("param1", particle_param1);
pp.query("param2", particle_param2);
pp.query("param3", particle_param3);
pp.query("verbose", particle_verbose);
}

} // namespace quokka
Expand Down
67 changes: 58 additions & 9 deletions src/problems/Gravity3D/gravity_3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,18 @@
struct BinaryOrbit {
};

// This is an ad-hoc test of particle creation and destruction.
// The initial condition consists of 2 particles with a mass of 1.0.
// In the first time step, 2^3 * 2 particles are created. Half of them are low-mass particles
// marked for destruction. In the next time step, low-mass particles are destroyed. There
// are 8 of them. Finally, in the third and last step, 2^3 * 2 particles are created.
// The final number of particles is 26

constexpr int particle_per_cell = 2;
constexpr int particle_spacing = 30;
constexpr double particle_low_mass = 1.0e-20; // very low mass particles marked for destruction
constexpr double dt_ = 0.001;
constexpr int n_particle_last = 26; // initial 2, 2^3 * 2 * 2 created, 8 destroyed

template <> struct quokka::EOS_Traits<BinaryOrbit> {
static constexpr double gamma = 1.0; // isothermal
Expand Down Expand Up @@ -66,7 +77,7 @@ template <> struct Physics_Traits<BinaryOrbit> {

namespace quokka
{
// Specialization for CIC particles
// Specialization for CIC particle creation
template <> struct ParticleCreationTraits<ParticleType::CIC> {
// Specialized nested ParticleChecker for CIC particles
template <typename problem_t> struct ParticleChecker {
Expand All @@ -79,13 +90,14 @@ template <> struct ParticleCreationTraits<ParticleType::CIC> {
// A simple demonstration of particle creation
// Could check density threshold or other state-based conditions
amrex::ignore_unused(state_arr, dx);
const int spacing = 16;
const int spacing = particle_spacing;
const bool is_create_particle_1 = current_time <= param1 && current_time + dt > param1;
const bool is_create_particle_2 = current_time <= param2 && current_time + dt > param2;
return ((is_create_particle_1 || is_create_particle_2) && (i != 0 && i % spacing == 0) && (j != 0 && j % spacing == 0) &&
(k != 0 && k % spacing == 0))
? particle_per_cell
: 0;
if ((is_create_particle_1 || is_create_particle_2) && (i != 0 && i % spacing == 0) && (j != 0 && j % spacing == 0) &&
(k != 0 && k % spacing == 0)) {
return particle_per_cell;
}
return 0;
}
};

Expand Down Expand Up @@ -113,7 +125,12 @@ template <> struct ParticleCreationTraits<ParticleType::CIC> {
const amrex::Real cell_volume = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
const amrex::Real cell_density = state_arr(i, j, k, HydroSystem<problem_t>::density_index);
const amrex::Real cell_mass = cell_density * cell_volume;
const amrex::Real particle_mass = 0.5 * cell_mass / num_particles; // Divide mass among particles
amrex::Real particle_mass = 0.5 * cell_mass / num_particles; // Divide mass among particles

// mark half of the particles as low-mass particles which will be destroyed in the next time step
if (i <= particle_spacing) {
particle_mass = particle_low_mass;
}

const amrex::Real vx = state_arr(i, j, k, HydroSystem<problem_t>::x1Momentum_index) / cell_density;
const amrex::Real vy = state_arr(i, j, k, HydroSystem<problem_t>::x2Momentum_index) / cell_density;
Expand Down Expand Up @@ -155,6 +172,37 @@ template <> struct ParticleCreationTraits<ParticleType::CIC> {
current_time, dt);
}
};

// Specialization for CIC particles destruction
template <> struct ParticleDestructionTraits<ParticleType::CIC> {
// Default nested ParticleChecker - determines if a particle should be destroyed
template <typename problem_t> struct ParticleChecker {
amrex::Real t_destroy = particle_param3;

// AMREX_GPU_HOST_DEVICE ParticleChecker(amrex::Real param1) : param1(param1) {}

template <typename ParticleType>
AMREX_GPU_DEVICE auto operator()(ParticleType &p, int mass_idx, amrex::Real current_time, amrex::Real dt) const -> bool
{
// Default implementation: destroy particles with mass < 1.0
amrex::ignore_unused(current_time, dt);

const bool is_small_mass = (p.rdata(mass_idx) < 2.0 * particle_low_mass);
const bool is_time = (current_time <= t_destroy && current_time + dt > t_destroy);
return is_small_mass && is_time;
}
};

// Main method to destroy particles - uses the helper implementation
template <typename problem_t, typename ContainerType>
static void destroyParticles(ContainerType *container, int mass_idx, int lev, amrex::Real current_time, amrex::Real dt)
{
// Use the common implementation with our checker type
ParticleDestructionImpl::destroyParticlesImpl<problem_t, ContainerType, ParticleDestructionTraits<ParticleType::CIC>::template ParticleChecker>(
container, mass_idx, lev, current_time, dt);
}
};

} // namespace quokka

template <> void QuokkaSimulation<BinaryOrbit>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
Expand Down Expand Up @@ -215,7 +263,8 @@ auto problem_main() -> int
// Problem initialization
QuokkaSimulation<BinaryOrbit> sim(BCs_cc);
sim.doPoissonSolve_ = 1; // enable self-gravity
// sim.initDt_ = 1.0e-2; // s
sim.initDt_ = dt_;
sim.maxDt_ = dt_;

// initialize
sim.setInitialConditions();
Expand All @@ -232,7 +281,7 @@ auto problem_main() -> int
double position_error = 0.0;
double position_norm = 0.0;

const int n_particle_expect = 2 + (3 * 3 * 3) * 2 * particle_per_cell; // 2 particles from the initial condition, (3*3*3)*2 particles from the creator
const int n_particle_expect = n_particle_last;

int status = 0; // Initialize to success

Expand Down
Loading

0 comments on commit 306cf35

Please sign in to comment.