From 731bfcb821a2d98c83c783a79f264b3ad3b8e327 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 21:08:47 +1100 Subject: [PATCH 01/10] add particle_destruction; gravity_3d compiles --- src/particles/PhysicsParticles.hpp | 45 +++++++++++--- src/particles/particle_creation.hpp | 2 +- src/particles/particle_destruction.hpp | 84 ++++++++++++++++++++++++++ 3 files changed, 120 insertions(+), 11 deletions(-) create mode 100644 src/particles/particle_destruction.hpp diff --git a/src/particles/PhysicsParticles.hpp b/src/particles/PhysicsParticles.hpp index 44b5250ef..4bfffb7c7 100644 --- a/src/particles/PhysicsParticles.hpp +++ b/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> = 0; @@ -71,6 +75,7 @@ class PhysicsParticleDescriptorBase 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 +92,8 @@ template [[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) { } @@ -293,6 +298,14 @@ template current_time, dt); } + void destroyParticles(int lev, amrex::Real current_time, amrex::Real dt) override + { + if (container_ != nullptr) { + ParticleDestructionTraits::template destroyParticles( + 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 +435,22 @@ template class PhysicsParticleRegister // Register a new particle type with specified properties template 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 descriptor; // Create the appropriate descriptor based on the particle type if (type == ParticleType::Rad) { descriptor = std::make_unique>( - 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>( - 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>( - 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 +560,23 @@ template 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)); diff --git a/src/particles/particle_creation.hpp b/src/particles/particle_creation.hpp index 4bfcb787c..eb37d3067 100644 --- a/src/particles/particle_creation.hpp +++ b/src/particles/particle_creation.hpp @@ -14,7 +14,7 @@ template 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 particle_checker; diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp new file mode 100644 index 000000000..1e8b1095f --- /dev/null +++ b/src/particles/particle_destruction.hpp @@ -0,0 +1,84 @@ +#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 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 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 to actually remove the invalid particles + container->Redistribute(lev); + } + } +} +} // namespace ParticleDestructionImpl + +// Traits class for specializing particle destruction behavior +template struct ParticleDestructionTraits { + // Default nested ParticleChecker - determines if a particle should be destroyed + template struct ParticleChecker { + AMREX_GPU_HOST_DEVICE ParticleChecker() = default; + + template + 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); + + // Check if mass is below threshold + return (p.rdata(mass_idx) < 1.0); // Destroy this particle + } + }; + + // Main method to destroy particles - uses the helper implementation + template + 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::template ParticleChecker>( + container, mass_idx, lev, current_time, dt); + } +}; + +} // namespace quokka + +#endif // PARTICLE_DESTRUCTION_HPP_ From 65e269a63ced43d635a0dbcf4884fdee916c74a4 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 22:13:54 +1100 Subject: [PATCH 02/10] add printParticleStatistics --- src/particles/PhysicsParticles.hpp | 19 +++++++++++++++++++ src/particles/particle_types.hpp | 2 ++ src/simulation.hpp | 11 +++++++++++ 3 files changed, 32 insertions(+) diff --git a/src/particles/PhysicsParticles.hpp b/src/particles/PhysicsParticles.hpp index 4bfffb7c7..a7565d5f9 100644 --- a/src/particles/PhysicsParticles.hpp +++ b/src/particles/PhysicsParticles.hpp @@ -70,6 +70,7 @@ 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 &rhs, int finest_lev, amrex::Real Gconst) = 0; virtual void driftParticles(int lev, amrex::Real dt) const = 0; @@ -214,6 +215,15 @@ template 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(container_->TotalNumberOfParticles(true, false)); + } + return 0; + } + #if AMREX_SPACEDIM == 3 // Implementation of mass deposition from particles to grid @@ -587,6 +597,15 @@ template class PhysicsParticleRegister } #endif // AMREX_SPACEDIM == 3 + // Print particle statistics + void printParticleStatistics() const + { + for (const auto &[type, descriptor] : particleRegistry_) { + amrex::Print() << "Particle type: " << getParticleTypeName(type); + amrex::Print() << ", Number of particles: " << descriptor->getNumParticles() << "\n"; + } + } + // Prevent copying or moving of the registry to ensure single ownership PhysicsParticleRegister(const PhysicsParticleRegister &) = delete; auto operator=(const PhysicsParticleRegister &) -> PhysicsParticleRegister & = delete; diff --git a/src/particles/particle_types.hpp b/src/particles/particle_types.hpp index 94d217878..152deb4d5 100644 --- a/src/particles/particle_types.hpp +++ b/src/particles/particle_types.hpp @@ -78,6 +78,7 @@ enum class ParticleType { // their own copies. inline amrex::Real particle_param1 = -1.0; // NOLINT inline amrex::Real particle_param2 = -1.0; // NOLINT +inline int particle_verbose = 0; // NOLINT print particle logistics //-------------------- Radiation particles -------------------- @@ -166,6 +167,7 @@ inline void particleParmParse() const amrex::ParmParse pp("particles"); pp.query("param1", particle_param1); pp.query("param2", particle_param2); + pp.query("verbose", particle_verbose); } } // namespace quokka diff --git a/src/simulation.hpp b/src/simulation.hpp index debff2071..3d990796a 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -323,6 +323,7 @@ template class AMRSimulation : public amrex::AmrCore void createDiagnostics(); void updateDiagnostics(); void doDiagnostics(); + void printParticleStatistics(); void WriteMetadataFile(std::string const &MetadataFileName) const; void ReadMetadataFile(std::string const &chkfilename); void WriteStatisticsFile(); @@ -1046,6 +1047,11 @@ template void AMRSimulation::evolve() WriteProjectionPlotfile(); } + // print particle statistics + if (quokka::particle_verbose > 0) { + printParticleStatistics(); + } + // write diagnostics doDiagnostics(); @@ -2364,6 +2370,11 @@ template void AMRSimulation::doDiagnostics() } } +template void AMRSimulation::printParticleStatistics() +{ + particleRegister_.printParticleStatistics(); +} + // do in-situ rendering with Ascent #ifdef AMREX_USE_ASCENT template void AMRSimulation::AscentCustomActions(conduit::Node const &blueprintMesh) From bc8f153da813c76784fac37a91002d631e5dc3c1 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 22:15:35 +1100 Subject: [PATCH 03/10] implement a test for particle destruction; passing test --- src/particles/particle_destruction.hpp | 11 +++++++++-- src/problems/Gravity3D/gravity_3d.cpp | 26 ++++++++++++++++++-------- src/simulation.hpp | 5 ++++- tests/Gravity3D.in | 5 +++-- 4 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp index 1e8b1095f..3d30197e9 100644 --- a/src/particles/particle_destruction.hpp +++ b/src/particles/particle_destruction.hpp @@ -55,7 +55,9 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev template struct ParticleDestructionTraits { // Default nested ParticleChecker - determines if a particle should be destroyed template struct ParticleChecker { - AMREX_GPU_HOST_DEVICE ParticleChecker() = default; + // amrex::Real param1 = ; + + // AMREX_GPU_HOST_DEVICE ParticleChecker(amrex::Real param1) : param1(param1) {} template AMREX_GPU_DEVICE auto operator()(ParticleType& p, int mass_idx, amrex::Real current_time, amrex::Real dt) const -> bool @@ -64,7 +66,12 @@ template struct ParticleDestructionTraits { amrex::ignore_unused(current_time, dt); // Check if mass is below threshold - return (p.rdata(mass_idx) < 1.0); // Destroy this particle + // return (p.rdata(mass_idx) < 1.0); // Destroy this particle + + const double t_destroy = 0.0015; + const bool is_small_mass = (p.rdata(mass_idx) < 2.0e-20); + const bool is_time = (current_time <= t_destroy && current_time + dt > t_destroy); + return is_small_mass && is_time; } }; diff --git a/src/problems/Gravity3D/gravity_3d.cpp b/src/problems/Gravity3D/gravity_3d.cpp index 24a649fc3..ecbd20d96 100644 --- a/src/problems/Gravity3D/gravity_3d.cpp +++ b/src/problems/Gravity3D/gravity_3d.cpp @@ -26,6 +26,9 @@ struct BinaryOrbit { }; 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; template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.0; // isothermal @@ -79,13 +82,14 @@ template <> struct ParticleCreationTraits { // 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 +117,12 @@ template <> struct ParticleCreationTraits { 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::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::x1Momentum_index) / cell_density; const amrex::Real vy = state_arr(i, j, k, HydroSystem::x2Momentum_index) / cell_density; @@ -215,7 +224,8 @@ auto problem_main() -> int // Problem initialization QuokkaSimulation 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 +242,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 = 2 + (2 * 2 * 2) * 2 * particle_per_cell; // 2 particles from the initial condition, (3*3*3)*2 particles from the creator int status = 0; // Initialize to success diff --git a/src/simulation.hpp b/src/simulation.hpp index 3d990796a..c0cc7118f 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -1013,6 +1013,9 @@ template void AMRSimulation::evolve() // Use the new type-aware particle creation method particleRegister_.createParticlesFromState(state_new_cc_[0], 0, cur_time, dt_[0]); + + // Use the new type-aware particle destruction method + particleRegister_.destroyParticles(0, cur_time, dt_[0]); #endif cur_time += dt_[0]; @@ -2154,7 +2157,7 @@ template void AMRSimulation::InitPhyParticles() CICParticles->SetVerbose(0); // Register with particle register - CIC particles allow creation - particleRegister_.registerParticleType(quokka::ParticleType::CIC, quokka::CICParticleMassIdx, -1, -1, false, true, CICParticles.get()); + particleRegister_.registerParticleType(quokka::ParticleType::CIC, quokka::CICParticleMassIdx, -1, -1, false, true, CICParticles.get(), true); // Initialize particles through derived class createInitialCICParticles(); diff --git a/tests/Gravity3D.in b/tests/Gravity3D.in index a6d0f8b08..d9eedaf1b 100644 --- a/tests/Gravity3D.in +++ b/tests/Gravity3D.in @@ -28,7 +28,8 @@ stop_time = 5.0 plotfile_interval = -1 checkpoint_interval = -1 -particles.param1 = 0.007 # create particles at this time, which is step 2 -particles.param2 = 0.012 # create particles at this time, which is step 3 +particles.param1 = 0.0005 # create particles at this time, which is step 0 +particles.param2 = 0.0025 # create particles at this time, which is step 2 +particles.verbose = 1 max_timesteps = 3 From 8dec836818fad6ef7374bf654884bad513107881 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 22:23:32 +1100 Subject: [PATCH 04/10] cleanup particle_destruction.hpp --- src/particles/particle_destruction.hpp | 15 ++-------- src/problems/Gravity3D/gravity_3d.cpp | 41 ++++++++++++++++++++++++-- 2 files changed, 42 insertions(+), 14 deletions(-) diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp index 3d30197e9..3cad0b7d4 100644 --- a/src/particles/particle_destruction.hpp +++ b/src/particles/particle_destruction.hpp @@ -55,24 +55,15 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev template struct ParticleDestructionTraits { // Default nested ParticleChecker - determines if a particle should be destroyed template struct ParticleChecker { - // amrex::Real param1 = ; - - // AMREX_GPU_HOST_DEVICE ParticleChecker(amrex::Real param1) : param1(param1) {} template 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); - - // Check if mass is below threshold - // return (p.rdata(mass_idx) < 1.0); // Destroy this particle - - const double t_destroy = 0.0015; - const bool is_small_mass = (p.rdata(mass_idx) < 2.0e-20); - const bool is_time = (current_time <= t_destroy && current_time + dt > t_destroy); - return is_small_mass && is_time; + amrex::ignore_unused(p, mass_idx, current_time, dt); + return false; } + }; // Main method to destroy particles - uses the helper implementation diff --git a/src/problems/Gravity3D/gravity_3d.cpp b/src/problems/Gravity3D/gravity_3d.cpp index ecbd20d96..0c18013ee 100644 --- a/src/problems/Gravity3D/gravity_3d.cpp +++ b/src/problems/Gravity3D/gravity_3d.cpp @@ -29,6 +29,7 @@ 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 { static constexpr double gamma = 1.0; // isothermal @@ -69,7 +70,7 @@ template <> struct Physics_Traits { namespace quokka { -// Specialization for CIC particles +// Specialization for CIC particle creation template <> struct ParticleCreationTraits { // Specialized nested ParticleChecker for CIC particles template struct ParticleChecker { @@ -164,6 +165,42 @@ template <> struct ParticleCreationTraits { current_time, dt); } }; + +// Specialization for CIC particles destruction +template <> struct ParticleDestructionTraits { + // Default nested ParticleChecker - determines if a particle should be destroyed + template struct ParticleChecker { + // amrex::Real param1 = ; + + // AMREX_GPU_HOST_DEVICE ParticleChecker(amrex::Real param1) : param1(param1) {} + + template + 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); + + // Check if mass is below threshold + // return (p.rdata(mass_idx) < 1.0); // Destroy this particle + + const double t_destroy = 0.0015; + const bool is_small_mass = (p.rdata(mass_idx) < 2.0e-20); + 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 + 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::template ParticleChecker>( + container, mass_idx, lev, current_time, dt); + } +}; + } // namespace quokka template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) @@ -242,7 +279,7 @@ auto problem_main() -> int double position_error = 0.0; double position_norm = 0.0; - const int n_particle_expect = 2 + (2 * 2 * 2) * 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 From 0b12e912bc8172e8c29e8474c322beb3f34b0bd7 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 22:36:24 +1100 Subject: [PATCH 05/10] add particle_param3 --- src/particles/particle_types.hpp | 2 ++ src/problems/Gravity3D/gravity_3d.cpp | 26 ++++++++++++++------------ tests/Gravity3D.in | 1 + 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/particles/particle_types.hpp b/src/particles/particle_types.hpp index 152deb4d5..5bfe0dd85 100644 --- a/src/particles/particle_types.hpp +++ b/src/particles/particle_types.hpp @@ -78,6 +78,7 @@ 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 -------------------- @@ -167,6 +168,7 @@ 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); } diff --git a/src/problems/Gravity3D/gravity_3d.cpp b/src/problems/Gravity3D/gravity_3d.cpp index 0c18013ee..7a76abd25 100644 --- a/src/problems/Gravity3D/gravity_3d.cpp +++ b/src/problems/Gravity3D/gravity_3d.cpp @@ -25,6 +25,13 @@ 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 @@ -87,7 +94,7 @@ template <> struct ParticleCreationTraits { 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; if ((is_create_particle_1 || is_create_particle_2) && (i != 0 && i % spacing == 0) && (j != 0 && j % spacing == 0) && - (k != 0 && k % spacing == 0)) { + (k != 0 && k % spacing == 0)) { return particle_per_cell; } return 0; @@ -119,7 +126,7 @@ template <> struct ParticleCreationTraits { const amrex::Real cell_density = state_arr(i, j, k, HydroSystem::density_index); const amrex::Real cell_mass = cell_density * cell_volume; 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; @@ -170,21 +177,17 @@ template <> struct ParticleCreationTraits { template <> struct ParticleDestructionTraits { // Default nested ParticleChecker - determines if a particle should be destroyed template struct ParticleChecker { - // amrex::Real param1 = ; + amrex::Real t_destroy = particle_param3; // AMREX_GPU_HOST_DEVICE ParticleChecker(amrex::Real param1) : param1(param1) {} template - AMREX_GPU_DEVICE auto operator()(ParticleType& p, int mass_idx, amrex::Real current_time, amrex::Real dt) const -> bool + 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); - - // Check if mass is below threshold - // return (p.rdata(mass_idx) < 1.0); // Destroy this particle - const double t_destroy = 0.0015; - const bool is_small_mass = (p.rdata(mass_idx) < 2.0e-20); + 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; } @@ -195,9 +198,8 @@ template <> struct ParticleDestructionTraits { 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::template ParticleChecker>( - container, mass_idx, lev, current_time, dt); + ParticleDestructionImpl::destroyParticlesImpl::template ParticleChecker>( + container, mass_idx, lev, current_time, dt); } }; diff --git a/tests/Gravity3D.in b/tests/Gravity3D.in index d9eedaf1b..448caebb9 100644 --- a/tests/Gravity3D.in +++ b/tests/Gravity3D.in @@ -30,6 +30,7 @@ checkpoint_interval = -1 particles.param1 = 0.0005 # create particles at this time, which is step 0 particles.param2 = 0.0025 # create particles at this time, which is step 2 +particles.param3 = 0.0015 # destroy particles at this time, which is step 1 particles.verbose = 1 max_timesteps = 3 From 25460062b42350a032d258660b19da505dfc4942 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Mar 2025 11:40:59 +0000 Subject: [PATCH 06/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/particles/PhysicsParticles.hpp | 9 +++++---- src/particles/particle_destruction.hpp | 24 +++++++++++------------- src/particles/particle_types.hpp | 2 +- src/simulation.hpp | 5 +---- 4 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/particles/PhysicsParticles.hpp b/src/particles/PhysicsParticles.hpp index a7565d5f9..30d4a8974 100644 --- a/src/particles/PhysicsParticles.hpp +++ b/src/particles/PhysicsParticles.hpp @@ -33,7 +33,7 @@ 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 + 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, bool allows_destruction = false) @@ -93,7 +93,8 @@ template [[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, bool allows_destruction = false) + 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) { } @@ -311,8 +312,8 @@ template void destroyParticles(int lev, amrex::Real current_time, amrex::Real dt) override { if (container_ != nullptr) { - ParticleDestructionTraits::template destroyParticles( - container_, this->getMassIndex(), lev, current_time, dt); + ParticleDestructionTraits::template destroyParticles(container_, this->getMassIndex(), lev, + current_time, dt); } } diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp index 3cad0b7d4..ac5aba00b 100644 --- a/src/particles/particle_destruction.hpp +++ b/src/particles/particle_destruction.hpp @@ -20,20 +20,20 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev // Iterate through all particles at this level for (typename ContainerType::ParIterType pti(*container, lev); pti.isValid(); ++pti) { - auto& particles = pti.GetArrayOfStructs(); + 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(); - + auto *parray = particles().data(); + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { - auto& p = parray[i]; // NOLINT - + 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 @@ -43,7 +43,7 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev } }); } - + // Redistribute particles to actually remove the invalid particles container->Redistribute(lev); } @@ -57,13 +57,12 @@ template struct ParticleDestructionTraits { template struct ParticleChecker { template - AMREX_GPU_DEVICE auto operator()(ParticleType& p, int mass_idx, amrex::Real current_time, amrex::Real dt) const -> bool + 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 @@ -71,9 +70,8 @@ template struct ParticleDestructionTraits { 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::template ParticleChecker>( - container, mass_idx, lev, current_time, dt); + ParticleDestructionImpl::destroyParticlesImpl::template ParticleChecker>( + container, mass_idx, lev, current_time, dt); } }; diff --git a/src/particles/particle_types.hpp b/src/particles/particle_types.hpp index 5bfe0dd85..c506350b8 100644 --- a/src/particles/particle_types.hpp +++ b/src/particles/particle_types.hpp @@ -79,7 +79,7 @@ enum class ParticleType { 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 +inline int particle_verbose = 0; // NOLINT print particle logistics //-------------------- Radiation particles -------------------- diff --git a/src/simulation.hpp b/src/simulation.hpp index c0cc7118f..2589816e7 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -2373,10 +2373,7 @@ template void AMRSimulation::doDiagnostics() } } -template void AMRSimulation::printParticleStatistics() -{ - particleRegister_.printParticleStatistics(); -} +template void AMRSimulation::printParticleStatistics() { particleRegister_.printParticleStatistics(); } // do in-situ rendering with Ascent #ifdef AMREX_USE_ASCENT From db00d676e4335e6953a37032f454815b186178a0 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 6 Mar 2025 22:51:14 +1100 Subject: [PATCH 07/10] mod printParticleStatistics --- src/particles/PhysicsParticles.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/particles/PhysicsParticles.hpp b/src/particles/PhysicsParticles.hpp index a7565d5f9..9b572ea2b 100644 --- a/src/particles/PhysicsParticles.hpp +++ b/src/particles/PhysicsParticles.hpp @@ -600,9 +600,10 @@ template class PhysicsParticleRegister // 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() << "Particle type: " << getParticleTypeName(type); - amrex::Print() << ", Number of particles: " << descriptor->getNumParticles() << "\n"; + amrex::Print() << getParticleTypeName(type) << ", " << descriptor->getNumParticles() << "\n"; } } From c6cf1fb25596d71ba37106b7f39146ac5de70389 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 7 Mar 2025 10:29:05 +1100 Subject: [PATCH 08/10] add TODO notes --- src/simulation.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/simulation.hpp b/src/simulation.hpp index 2589816e7..950372809 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -1012,9 +1012,11 @@ template void AMRSimulation::evolve() kickParticlesAllLevels(dt_[0]); // Use the new type-aware particle creation method + // TODO(cch): Need to take care of AMR subscycling particleRegister_.createParticlesFromState(state_new_cc_[0], 0, cur_time, dt_[0]); // Use the new type-aware particle destruction method + // TODO(cch): Need to take care of AMR subscycling particleRegister_.destroyParticles(0, cur_time, dt_[0]); #endif From ad73565e5e2564e2ea523fb385ca4a4f7771be9b Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 7 Mar 2025 11:31:51 +1100 Subject: [PATCH 09/10] comments about amr subcycling and particle destruction --- src/particles/particle_destruction.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp index ac5aba00b..9fa4389a8 100644 --- a/src/particles/particle_destruction.hpp +++ b/src/particles/particle_destruction.hpp @@ -44,7 +44,9 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev }); } - // Redistribute particles to actually remove the invalid particles + // 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); } } From e9d2860e5dea7e042b92758040dabf1d24f32bf3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 7 Mar 2025 00:32:18 +0000 Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/particles/particle_destruction.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/particles/particle_destruction.hpp b/src/particles/particle_destruction.hpp index 9fa4389a8..5d7eecde0 100644 --- a/src/particles/particle_destruction.hpp +++ b/src/particles/particle_destruction.hpp @@ -45,8 +45,8 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev } // 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. + // 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); } }