Add print of SN explosion counts (quokka-astro#1546)

### Description
Print the number of stars that went SN explosion in each step (if >0)
when `verbose = true`.

Also change the format of some print logs in the particle module to make
them more consistent and avoid excessive logs.

### Related issues
Are there any GitHub issues that are fixed by this pull request? Add a
link to them here.

### 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:_
- [ ] I have added a description (see above).
- [ ] I have added a link to any related issues (if applicable; see
above).
- [ ] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [ ] 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>

Author: chongchonghe

src/particles/PhysicsParticles.hpp

@@ -145,9 +145,9 @@ class PhysicsParticleDescriptorBase
 	//----- Methods that are implemented for some but not all particle types, so they cannot be pure virtual -----
 
 	virtual auto depositSN(amrex::MultiFab & /*state*/, std::array<amrex::MultiFab, AMREX_SPACEDIM> const * /*state_fc*/, int /*lev*/, amrex::Real /*time*/,
-			       amrex::Real /*dt*/) -> amrex::Real
+			       amrex::Real /*dt*/) -> std::pair<int, amrex::Real>
 	{
-		return 0.0_rt;
+		return {0, 0.0_rt};
 	}
 
 	virtual void computeSinkAccretion(amrex::MultiFab &state, amrex::MultiFab &state_accretion_rate,
@@ -610,8 +610,9 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
 
 	// Implementation of supernova energy and momentum deposition from particles to grid
 	auto depositSN(amrex::MultiFab &state, std::array<amrex::MultiFab, AMREX_SPACEDIM> const *state_fc, int lev, amrex::Real time, amrex::Real dt)
-	    -> amrex::Real override
+	    -> std::pair<int, amrex::Real> override
 	{
+		int num_sn_explosions = 0;
 		amrex::Real max_velocity = 0.0;
 
 		if (this->container_ != nullptr && this->getEvolutionStageIndex() >= 0) {
@@ -621,16 +622,19 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
 								 "UnitSystem must be CGS for particleMeshInteraction");
 
 				// Deposit supernova energy and momentum from all particles. This also updates the evolution stage of the particles.
-				max_velocity = SNDeposition<ContainerType, problem_t>(this->container_, state, state_fc, lev, time, dt, this->getMassIndex(),
-										      this->getEvolutionStageIndex(), this->getBirthTimeIndex());
+				auto [sn_count, vel] =
+				    SNDeposition<ContainerType, problem_t>(this->container_, state, state_fc, lev, time, dt, this->getMassIndex(),
+									   this->getEvolutionStageIndex(), this->getBirthTimeIndex());
+				num_sn_explosions = sn_count;
+				max_velocity = vel;
 			} else {
 				// Only update evolution stage but not deposit energy/momentum
 				SNFeedbackUtils::updateEvolutionStage(this->container_, lev, time + dt, this->getBirthTimeIndex(),
 								      this->getEvolutionStageIndex());
 			}
 		}
 
-		return max_velocity;
+		return {num_sn_explosions, max_velocity};
 	}
 
 	// compute accretion rate
@@ -807,16 +811,18 @@ template <typename problem_t> class PhysicsParticleRegister
 
 	// Deposit supernova energy and momentum from all particles
 	auto depositSN(amrex::MultiFab &state, std::array<amrex::MultiFab, AMREX_SPACEDIM> const *state_fc, int lev, amrex::Real time, amrex::Real dt)
-	    -> amrex::Real
+	    -> std::pair<int, amrex::Real>
 	{
 		const BL_PROFILE("PhysicsParticleRegister::depositSN()");
+		int total_sn_explosions = 0;
 		amrex::Real max_velocity = 0.0;
 		// Each particle type handles its own buffer creation and roundoff independently
 		for (const auto &[type, descriptor] : particleRegistry_) {
-			const amrex::Real max_velocity_ = descriptor->depositSN(state, state_fc, lev, time, dt);
-			max_velocity = std::max(max_velocity, max_velocity_);
+			auto [sn_count, velocity] = descriptor->depositSN(state, state_fc, lev, time, dt);
+			total_sn_explosions += sn_count;
+			max_velocity = std::max(max_velocity, velocity);
 		}
-		return max_velocity;
+		return {total_sn_explosions, max_velocity};
 	}
 
 	// Implementation of computeSinkAccretion
@@ -898,7 +904,7 @@ template <typename problem_t> class PhysicsParticleRegister
 				}
 			}
 			if (!found) {
-				amrex::Print() << "Warning: Requested particle type '" << requestedName << "' is not registered.\n";
+				amrex::Print() << "[WARNING] Requested particle type '" << requestedName << "' is not registered.\n";
 			}
 		}
 	}
@@ -1020,7 +1026,7 @@ template <typename problem_t> class PhysicsParticleRegister
 	void printParticleStatistics() const
 	{
 		const BL_PROFILE("PhysicsParticleRegister::printParticleStatistics()");
-		amrex::Print() << ">>> Particle statistics:\n";
+		amrex::Print() << "[PARTICLES] Statistics:\n";
 		amrex::Print() << fmt::format("{:<20}{:>15}\n", "Particle type", "Number of particles");
 
 		for (const auto &[type, descriptor] : particleRegistry_) {

src/particles/particle_deposition.hpp

@@ -272,7 +272,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void depositThermalKineticMomentumSNR(
 
 template <typename ContainerType, typename problem_t>
 void depositToBuffer(ContainerType *container, amrex::MultiFab &state, amrex::MultiFab &state_buffer, int lev, amrex::Real time, amrex::Real dt, int mass_index,
-		     int evolutionStageIndex, int birthTimeIndex, const SNScheme SN_scheme_d)
+		     int evolutionStageIndex, int birthTimeIndex, const SNScheme SN_scheme_d, int *p_sn_count = nullptr)
 {
 	const BL_PROFILE("SNFeedbackUtils::depositToBuffer()");
 	constexpr amrex::Real stencil_volume = 4.0 / 3.0 * M_PI * SN_stencil_size * SN_stencil_size * SN_stencil_size;
@@ -330,6 +330,11 @@ void depositToBuffer(ContainerType *container, amrex::MultiFab &state, amrex::Mu
 			const bool is_sn_progenitor = (p.idata(evolutionStageIndex) == static_cast<int>(StellarEvolutionStage::SNProgenitor));
 
 			if (is_sn_progenitor && step_end_time > p.rdata(birthTimeIndex + 1)) {
+				// Count this SN explosion
+				if (p_sn_count != nullptr) {
+					amrex::Gpu::Atomic::AddNoRet(p_sn_count, 1);
+				}
+
 				// Update the particle's evolution stage to SNRemnant
 				p.idata(evolutionStageIndex) = static_cast<int>(StellarEvolutionStage::SNRemnant);
 
@@ -621,7 +626,7 @@ void updateEvolutionStage(ContainerType *container, int lev_min, amrex::Real ste
 
 template <typename ContainerType, typename problem_t>
 auto SNDeposition(ContainerType *container, amrex::MultiFab &state, std::array<amrex::MultiFab, AMREX_SPACEDIM> const *state_fc, int lev, amrex::Real time,
-		  amrex::Real dt, int mass_index, int evolutionStageIndex, int birthTimeIndex) -> Real
+		  amrex::Real dt, int mass_index, int evolutionStageIndex, int birthTimeIndex) -> std::pair<int, Real>
 {
 	const BL_PROFILE("[particle_deposition] SNDeposition()");
 	static_assert(SN_stencil_size <= 3,
@@ -638,9 +643,13 @@ auto SNDeposition(ContainerType *container, amrex::MultiFab &state, std::array<a
 	amrex::Gpu::Buffer<amrex::Real> max_velocity_buffer({0.0});
 	amrex::Real *p_max_velocity = max_velocity_buffer.data();
 
+	// Initialize SN explosion count tracking
+	amrex::Gpu::Buffer<int> sn_count_buffer({0});
+	int *p_sn_count = sn_count_buffer.data();
+
 	// Step 1: Local deposition within each box
 	SNFeedbackUtils::depositToBuffer<ContainerType, problem_t>(container, state, state_buffer, lev, time, dt, mass_index, evolutionStageIndex,
-								   birthTimeIndex, SN_scheme_d);
+								   birthTimeIndex, SN_scheme_d, p_sn_count);
 
 	// Step 2: Sum boundary values
 	state_buffer.SumBoundary(container->Geom(lev).periodicity());
@@ -651,12 +660,16 @@ auto SNDeposition(ContainerType *container, amrex::MultiFab &state, std::array<a
 	// Step 3: Add the buffer to the state
 	SNFeedbackUtils::addBufferToState<problem_t>(state, state_fc, state_buffer, SN_scheme_d, p_max_velocity);
 
-	// Step 4: Check maximum velocity and print warning if needed
+	// Step 4: Check maximum velocity and SN count, then reduce across ranks
 	auto *h_max_velocity = max_velocity_buffer.copyToHost();
 	Real max_velocity = h_max_velocity[0];
 	amrex::ParallelDescriptor::ReduceRealMax(max_velocity);
 
-	return max_velocity;
+	auto *h_sn_count = sn_count_buffer.copyToHost();
+	int sn_count = h_sn_count[0];
+	amrex::ParallelDescriptor::ReduceIntSum(sn_count);
+
+	return {sn_count, max_velocity};
 }
 
 #endif // AMREX_SPACEDIM == 3

src/particles/particle_destruction.hpp

@@ -3,6 +3,7 @@
 
 #include "AMReX_BLProfiler.H"
 #include "particle_types.hpp"
+#include <fmt/format.h>
 
 namespace quokka
 {
@@ -79,8 +80,11 @@ static void destroyParticlesImpl(ContainerType *container, int mass_idx, int lev
 
 			// Print the total number of particles destroyed at this time step
 			if (amrex::ParallelDescriptor::IOProcessor()) {
-				amrex::Print() << ">>>Particle destruction:\n\tTime: " << current_time << " - Destroyed " << global_total_destroyed
-					       << " particles (from level " << lev_min << " and above)\n\n";
+				if (particle_verbose > 0) {
+					amrex::Print()
+					    << fmt::format("[PARTICLES] Particle destruction: Time: {} - Destroyed {} particles (from level {} and above)\n",
+							   current_time, global_total_destroyed, lev_min);
+				}
 			}
 		}
 	}

src/simulation.hpp

@@ -1922,12 +1922,17 @@ template <typename problem_t> void AMRSimulation<problem_t>::particleMeshInterac
 	particleRegister_.createParticlesFromState(state_new_cc_[lev], accretion_rate_at_level, lev, time, dt, state_fc_ptr, verbose);
 
 	// Deposit the SN particles into the MultiFab
-	const amrex::Real max_velocity = particleRegister_.depositSN(state_new_cc_[lev], state_fc_ptr, lev, time, dt);
+	const auto [num_sn_explosions, max_velocity] = particleRegister_.depositSN(state_new_cc_[lev], state_fc_ptr, lev, time, dt);
+
+	// Print SN explosion count if verbose and non-zero
+	if (verbose && num_sn_explosions > 0) {
+		amrex::Print() << fmt::format("[PARTICLES] SN explosions: Time: {} - {} stars went supernova at level {}\n", time, num_sn_explosions, lev);
+	}
 
 	// Check if the maximum velocity is greater than the threshold
 	constexpr amrex::Real v_over_c_threshold = 0.03;
 	if (max_velocity > v_over_c_threshold * C::c_light) {
-		amrex::Print() << "WARNING: SN remnant net velocity (" << max_velocity / C::c_light << " c) greater than " << v_over_c_threshold
+		amrex::Print() << "[WARNING] SN remnant net velocity (" << max_velocity / C::c_light << " c) greater than " << v_over_c_threshold
 			       << " c threshold!" << "\n";
 	}
 }