A new SN feedback model (#995)

### Description

In this PR I implement a new SNe feedback model that smoothly
transitions between thermal feedback, kinetic feedback, and momentum
feedback, depending on how well the cooling/shell formation phase of the
SNR evolution is resolved. An innovation of this scheme is that we do a
two stage deposition scheme to account for SNR overlapping and get rid
of depedence on the order that SNe events occur.

I have tested three models:
1. Model 0: Pure thermal deposition (for comparison only, not for
production)
2. Model 1: Thermal energy deposition in the resolved regime, thermal +
kinetic in the intermediate regime, and terminal momentum in the
unresolved regime.
3. Model 2: Kinetic energy deposition in the resolved regime, thermal +
kinetic in the intermediate regime, and terminal momentum in the
unresolved regime.

At a SN event, we distribute mass, momentum, and energy uniformly among
the cells that intercect with s sphere of radius `3 dx` centered at the
center of the cell that the SN is located in. No background smoothing is
conducted before deposition.

The results looks great! See figure below. Both model 1 and 2 produces
'correct' and almost identical results. The only difference is the
maximum temperature (defined as the maximum temperature in any cell
throughout the simulation) in model 2 is slightly lower than that in
model 1. In model 2, right after SN explosion, the gas is cool but moves
fast, and the shock wave immediately heat the gas to comparable
temperature to that in the thermal model. After some experiments, I am
convinced that the thermal model and kinetic model (in the resolved
regime) has no difference in accuracy or performance related to
timestep. This is because, imagine a star explodes into vacuum, q simple
calculation will show that the sound speed of the hot SNR in the pure
thermal model is about $1/\sqrt{3}$ times the bulk velocity of the cool
SNR in the pure kinetic model (microscopic thermal motion has three
degrees of freedom while kinetic motion has one). The thermal model
turns out to have slightly longer timestep overall. In the TIGRESS model
they advocate thermal model because they have over cooling issue, but we
don't have this issue at all (thanks to Ben's awesome Grackle cooling
module which has a robust ODE solver).


![compare](https://github.com/user-attachments/assets/23895e3d-e00a-46a5-af30-a62d70509a02)

> Reproducing Figure 6 of Kim & Ostriker 2017. The spacial resolution
and SN ejecta mass/energy/momentum in our simulations are the same as
theirs.

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

Author: chongchonghe

src/particles/PhysicsParticles.hpp

@@ -9,12 +9,12 @@
 #include <string>
 
 #include "AMReX_Array4.H"
+#include "AMReX_BLassert.H"
 #include "AMReX_MultiFab.H"
 #include "AMReX_ParticleInterpolators.H"
 #include "AMReX_REAL.H"
 #include "AMReX_SPACE.H"
 #include "AMReX_Vector.H"
-#include "hydro/hydro_system.hpp"
 #include "particle_creation.hpp"
 #include "particle_deposition.hpp"
 #include "particle_destruction.hpp"
@@ -118,7 +118,8 @@ class PhysicsParticleDescriptorBase
 	[[nodiscard]] virtual auto computeMaxParticleSpeed(int lev) const -> amrex::ValLocPair<amrex::Real, amrex::RealVect> = 0;
 
 	// Methods that are implemented for some but not all particle types, so they cannot be pure virtual
-	virtual void depositSN(const amrex::Vector<amrex::MultiFab *> &state, int lev_min, amrex::Real step_end_time) { /* Default empty implementation */ }
+	virtual void depositSN(amrex::MultiFab &state, amrex::MultiFab &state_buffer, int lev, amrex::Real time, amrex::Real dt)
+	{ /* Default empty implementation */ }
 #endif // AMREX_SPACEDIM == 3
 };
 
@@ -570,7 +571,7 @@ template <typename ContainerType, typename problem_t, ParticleType particleType>
 				}
 
 				// Get the units data for this particle type
-				const auto &unitsData = quokka::get_units_data();
+				const auto &unitsData = get_units_data();
 				if (unitsData.find(particleType_) == unitsData.end()) {
 					amrex::Abort(
 					    "Error: Particle type not defined in units data map. Please add units for this particle type in get_units_data().");
@@ -640,19 +641,21 @@ class StarParticleDescriptor : public PhysicsParticleDescriptor<ContainerType, p
 
 #if AMREX_SPACEDIM == 3
 	// Implementation of supernova energy and momentum deposition from particles to grid
-	void depositSN(const amrex::Vector<amrex::MultiFab *> &state, int lev_min, amrex::Real step_end_time) override
+	void depositSN(amrex::MultiFab &state, amrex::MultiFab &state_buffer, int lev, amrex::Real time, amrex::Real dt) override
 	{
 		if (this->container_ != nullptr && this->getEvolutionStageIndex() >= 0) {
-			// Deposit supernova energy and momentum from level lev_min to finest level
-			// zero_out_input is false because we want to accumulate supernova contributions
-			// vol_weight is false because SNDeposition does the volume weighting
-			amrex::ParticleToMesh(*this->container_, state, lev_min, -1,
-					      SNDeposition{step_end_time, this->getMassIndex(), HydroSystem<problem_t>::density_index,
-							   this->getBirthTimeIndex(), this->getEvolutionStageIndex()},
-					      false, false);
-
-			// Update particle evolution stages after deposition
-			updateEvolutionStage(this->container_, lev_min, step_end_time, this->getBirthTimeIndex(), this->getEvolutionStageIndex());
+			if (!quokka::disable_SN_feedback) {
+				// Requires CGS units
+				AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Physics_Traits<problem_t>::unit_system == UnitSystem::CGS,
+								 "UnitSystem must be CGS for particleMeshInteraction");
+
+				// Deposit supernova energy and momentum from all particles. This also updates the evolution stage of the particles.
+				SNDeposition<ContainerType, problem_t>(this->container_, state, state_buffer, lev, time, dt, this->getMassIndex(),
+								       this->getEvolutionStageIndex(), this->getBirthTimeIndex());
+			} else {
+				// Only update evolution stage but not deposit energy/momentum
+				updateEvolutionStage(this->container_, lev, time + dt, this->getBirthTimeIndex(), this->getEvolutionStageIndex());
+			}
 		}
 	}
 #endif // AMREX_SPACEDIM == 3
@@ -682,6 +685,17 @@ template <typename problem_t> class PhysicsParticleRegister
 		return false;
 	}
 
+	// Check if registry contains any star particles
+	[[nodiscard]] auto HasStarParticles() const -> bool
+	{
+		for (const auto &[name, descriptor] : particleRegistry_) {
+			if (descriptor->isStarParticle()) {
+				return true;
+			}
+		}
+		return false;
+	}
+
 	// Utility method to convert particle type to string name (for writing plotfiles/checkpoints)
 	[[nodiscard]] static auto getParticleTypeName(ParticleType type) -> std::string
 	{
@@ -786,13 +800,12 @@ template <typename problem_t> class PhysicsParticleRegister
 	}
 
 	// Deposit supernova energy and momentum from all particles
-	void depositSN(const amrex::Vector<amrex::MultiFab *> &state, int lev_min, amrex::Real step_end_time)
+	void depositSN(amrex::MultiFab &state, amrex::MultiFab &state_buffer, int lev, amrex::Real time, amrex::Real dt)
 	{
-		// this function is only implemented for particles that belong to the StarParticleDescriptor class whose unique signature is that the
-		// isStarParticle() method returns true
+		// this function is only implemented for some particle types, so we specify the particle type manually here
 		for (const auto &[type, descriptor] : particleRegistry_) {
 			if (descriptor->isStarParticle()) {
-				descriptor->depositSN(state, lev_min, step_end_time);
+				descriptor->depositSN(state, state_buffer, lev, time, dt);
 			}
 		}
 	}