Particle destructor

src/particles/PhysicsParticles.hpp

@@ -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"
 
@@ -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)
 	{
 	}
 
@@ -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;
@@ -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
 };
@@ -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)
 	{
 	}
 
@@ -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
@@ -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
 	{
@@ -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 {
@@ -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));
@@ -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;

src/particles/particle_creation.hpp

@@ -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;
 

src/particles/particle_destruction.hpp

@@ -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_

src/particles/particle_types.hpp

@@ -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 --------------------
 
@@ -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

src/problems/Gravity3D/gravity_3d.cpp

@@ -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
@@ -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 {
@@ -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;
 		}
 	};
 
@@ -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;
@@ -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)
@@ -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();
@@ -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