Modular stellar-evolution framework for Star particles

docs/markdown/particles.md

@@ -393,3 +393,41 @@ In each of the three tests, the simulation is run twice, one initially at rest a
 #### Random Blast Test
 
 The `RandomBlast` problem provides a testbed for multiple SN explosions with various ambient conditions and bulk flows.
+
+## Star Particle Type (modular stellar evolution)
+
+Star particles (`ParticleSwitch::Star`) represent individual stars whose radius and
+luminosity evolve through a **pluggable stellar-evolution model** selected at compile time.
+
+### Particle attributes
+
+Real components: `mass`, `vx`, `vy`, `vz`, `birth_time`, `mdot` (accretion rate, set by the
+accretion module), `radius` (set by the stellar model), and `lum` (luminosity per radiation
+group; occupies the last `nGroups` slots). A model may declare extra real/integer components
+via `nExtraReal` / `nExtraInt`.
+
+### Choosing a model
+
+The model is selected via `Particle_Traits<problem_t>::stellar_model`, which defaults to
+`ToyStellarModel` through `DefaultParticleTraits`. A problem can override it by adding
+`using stellar_model = MyModel;` to its `Particle_Traits` specialization. A model is a struct
+of GPU device functions; the dispatcher `StellarUpdate::updateStellarProperties` reads the
+particle's `mass` and `mdot`, calls the model, and stores `radius` and `lum` once per coarse
+step (operator-split, after accretion).
+
+### Toy model
+
+`ToyStellarModel` is stateless:
+
+- Radius: $R = R_\odot (M/M_\odot)^{0.4}$
+- Stellar luminosity: $L_\star = L_\odot (M/M_\odot)^{3.5}$
+- Accretion luminosity: $L_\mathrm{acc} = G M \dot{M} / R$
+- Stored luminosity: $L = L_\star + L_\mathrm{acc}$
+
+### Validation test
+
+`ParticleStarEvolution` places one Star particle in a uniform medium and lets it accrete via
+the grid Bondi accretion module (small-mass regime, $\dot{M}\approx$ const). It asserts, each
+step, that the particle's stored radius and luminosity match the closed-form laws above
+within ~2% (the tolerance absorbing the one-timestep lag between accretion and the stellar
+update), and reports the mean accretion rate against the analytic Bondi value.

inputs/ParticleStarEvolution.toml

@@ -0,0 +1,27 @@
+# Uniform-grid Bondi accretion onto a single Star particle (toy stellar evolution validation)
+# Domain sized so that r_B ~ 0.1 * dx, putting the accretion in the sub-grid Bondi regime.
+
+geometry.prob_lo = [-1.2e19, -1.2e19, -1.2e19]
+geometry.prob_hi = [ 1.2e19,  1.2e19,  1.2e19]
+quokka.bc = ["foextrap", "foextrap", "foextrap"]
+
+amr.n_cell = [64, 64, 64]
+amr.max_level = 0
+amr.blocking_factor = 16
+amr.max_grid_size = 32
+
+do_reflux = 1
+do_subcycle = 0
+do_tracers = 0
+
+max_timesteps = 10
+
+particles.verbose = 1
+
+problem.M0_in_Msun = 0.1
+problem.rho0 = 1.6726e-24  # ~ m_p (n_H ~ 1)
+problem.t_end_over_t_b = 300.0
+
+plotfile_interval = -1
+checkpoint_interval = -1
+tiny_profiler.enabled = 0

inputs/star.txt

@@ -0,0 +1,2 @@
+1
+0.0 0.0 0.0 1.988409871e+32 0.0 0.0 0.0 0.0

src/QuokkaSimulation.hpp

@@ -308,6 +308,7 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 	void createInitialCICRadParticles() override;
 	void createInitialStochasticStellarPopParticles() override;
 	void createInitialSinkParticles() override;
+	void createInitialStarParticles() override;
 	void createInitialTestParticles() override;
 #endif // AMREX_SPACEDIM == 3
 	void advanceSingleTimestepAtLevel(int lev, amrex::Real time, amrex::Real dt_lev, int ncycle) override;
@@ -954,6 +955,13 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::createInitialSin
 	// note: an implementation is only effective if Sink_particles are used
 }
 
+template <typename problem_t> void QuokkaSimulation<problem_t>::createInitialStarParticles()
+{
+	const BL_PROFILE("QuokkaSimulation::createInitialStarParticles()");
+	// Default: no Star particles. A problem overrides this to place particles at t=0.
+	// Only effective when ParticleSwitch::Star is enabled for the problem.
+}
+
 template <typename problem_t> void QuokkaSimulation<problem_t>::createInitialTestParticles()
 {
 	const BL_PROFILE("QuokkaSimulation::createInitialTestParticles()");

src/linear_advection/AdvectionSimulation.hpp

@@ -99,6 +99,7 @@ template <typename problem_t> class AdvectionSimulation : public AMRSimulation<p
 	void createInitialCICRadParticles() override;
 	void createInitialStochasticStellarPopParticles() override;
 	void createInitialSinkParticles() override;
+	void createInitialStarParticles() override;
 	void createInitialTestParticles() override;
 #endif // AMREX_SPACEDIM == 3
 	void advanceSingleTimestepAtLevel(int lev, amrex::Real time, amrex::Real dt_lev, int /*ncycle*/) override;
@@ -224,6 +225,12 @@ template <typename problem_t> void AdvectionSimulation<problem_t>::createInitial
 	// note: an implementation is only effective if Sink particles are used
 }
 
+template <typename problem_t> void AdvectionSimulation<problem_t>::createInitialStarParticles()
+{
+	// Optional implementation
+	// note: an implementation is only effective if Star particles are used
+}
+
 template <typename problem_t> void AdvectionSimulation<problem_t>::createInitialTestParticles()
 {
 	// Optional implementation

src/particles/PhysicsParticles.hpp

@@ -54,6 +54,8 @@ namespace quokka
 			return "StochasticStellarPop";
 		case ParticleType::Sink:
 			return "Sink";
+		case ParticleType::Star:
+			return "Star";
 		default:
 			AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown particle type");
 			return "Unknown";
@@ -75,6 +77,8 @@ namespace quokka
 			return "ParticleSwitch::StochasticStellarPop";
 		case ParticleType::Sink:
 			return "ParticleSwitch::Sink";
+		case ParticleType::Star:
+			return "ParticleSwitch::Star";
 		default:
 			AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown particle type");
 			return "ParticleSwitch::Unknown";
@@ -101,6 +105,9 @@ namespace quokka
 	if (particle_type_name == "Sink") {
 		return ParticleType::Sink;
 	}
+	if (particle_type_name == "Star") {
+		return ParticleType::Star;
+	}
 	return std::nullopt;
 }
 
@@ -887,6 +894,8 @@ template <typename problem_t> class PhysicsParticleRegister
 				return "StochasticStellarPop_particles";
 			case ParticleType::Sink:
 				return "Sink_particles";
+			case ParticleType::Star:
+				return "Star_particles";
 			default:
 				AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "Unknown particle type");
 				return "Unknown_particles";
@@ -918,6 +927,11 @@ template <typename problem_t> class PhysicsParticleRegister
 		} else if constexpr (particleType == ParticleType::Sink) {
 			descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::Sink>>(
 			    container, SinkParticleMassIdx, -1, -1, -1, true, false, -1, true, -1, SinkParticleMdotIdx, SinkParticleLxIdx);
+		} else if constexpr (particleType == ParticleType::Star) {
+			constexpr int lum_idx = (Physics_Traits<problem_t>::nGroups > 0) ? StarParticleLumIdx : -1;
+			descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::Star>>(
+			    container, StarParticleMassIdx, lum_idx, StarParticleBirthTimeIdx, -1, /*allows_creation=*/false,
+			    /*allows_destruction=*/false, /*evolution_stage_idx=*/-1, /*allows_accretion=*/true, /*mass_at_birth_idx=*/-1, StarParticleMdotIdx);
 		} else if constexpr (particleType == ParticleType::Test) {
 			descriptor = std::make_unique<PhysicsParticleDescriptor<ContainerType, problem_t, ParticleType::Test>>(
 			    container, TestParticleMassIdx, TestParticleLumIdx, TestParticleBirthTimeIdx, TestParticleDeathTimeIdx, true, true,

src/particles/particle_types.hpp

@@ -4,6 +4,7 @@
 #include "AMReX_AmrParticles.H"
 #include "AMReX_Enum.H"
 #include "AMReX_ParIter.H"
+#include "particles/stellar_models.hpp"
 #include "physics_info.hpp"
 
 // Function to create bit flags: bitflag(position) = 2^(position - 1)
@@ -22,7 +23,8 @@ enum class ParticleSwitch : unsigned int {
 	CICRad = bitflag<3>(),		     // Combined gravitating-radiating particles, = 0b0100
 	StochasticStellarPop = bitflag<4>(), // Stellar population particles, = 0b1000
 	Sink = bitflag<5>(),		     // Sink particles, = 0b10000
-	Test = bitflag<6>()		     // Test particles with all features enabled, = 0b100000
+	Test = bitflag<6>(),		     // Test particles with all features enabled, = 0b100000
+	Star = bitflag<7>()		     // Star particles, = 0b1000000
 };
 
 // Enable bitwise operations on the enum class
@@ -49,8 +51,14 @@ constexpr auto operator&(ParticleSwitch flags, ParticleSwitch flag) -> bool
 // };
 // - static constexpr TestEnum particle_switch = TestEnum::MISTAKE;
 // - static constexpr ParticleSwitch particle_switch = ParticleSwitch::CIC | TestEnum::MISTAKE;
-template <typename problem_t> struct Particle_Traits {
+struct DefaultParticleTraits {
 	static constexpr ParticleSwitch particle_switch = ParticleSwitch::None; // Determines which particle types are enabled using bitwise flags.
+	using stellar_model = quokka::ToyStellarModel;				// Default stellar-evolution model
+};
+
+// This struct should be specialized by the user application code to configure particle behavior.
+// Inherits from DefaultParticleTraits so that new parameters added to it don't require updating every problem.
+template <typename problem_t> struct Particle_Traits : DefaultParticleTraits {
 };
 
 // Static assertion helper to verify that particle_switch is of the correct type
@@ -75,7 +83,8 @@ enum class ParticleType {
 	CICRad,		      // Gravitating radiation particles
 	StochasticStellarPop, // Stellar population particles
 	Sink,		      // Sink particles
-	Test		      // Test particles with all features enabled
+	Test,		      // Test particles with all features enabled
+	Star		      // Star particles
 };
 
 // Compile-time particle metadata.
@@ -339,6 +348,20 @@ constexpr int SinkParticleRealComps = 8; // mass, vx, vy, vz, mdot, Lx, Ly, Lz
 using SinkParticleContainer = amrex::AmrParticleContainer<SinkParticleRealComps>;
 using SinkParticleIterator = amrex::ParIter<SinkParticleRealComps>;
 
+//-------------------- Star particles --------------------
+// The layout enum and index constants live in a standalone header so that stellar
+// models can include them without pulling in the full particle type machinery.
+
+#include "particles/star_particle_indices.H"
+
+// Number of components = fixed scalars + nGroups luminosity slots + model extras.
+template <typename problem_t>
+constexpr int StarParticleRealComps = StarParticleFixedComps + Physics_Traits<problem_t>::nGroups + Particle_Traits<problem_t>::stellar_model::nExtraReal;
+template <typename problem_t> constexpr int StarParticleIntComps = Particle_Traits<problem_t>::stellar_model::nExtraInt;
+
+template <typename problem_t> using StarParticleContainer = amrex::AmrParticleContainer<StarParticleRealComps<problem_t>, StarParticleIntComps<problem_t>>;
+template <typename problem_t> using StarParticleIterator = amrex::ParIter<StarParticleRealComps<problem_t>, StarParticleIntComps<problem_t>>;
+
 //-------------------- Component Names for I/O --------------------
 
 // Helper function to generate component names from an enum type
@@ -386,6 +409,8 @@ template <ParticleType particleType, typename problem_t> auto getParticleRealCom
 		return expandEnumNames<CICRadParticleRealIdx, CICRadParticleRealComps<problem_t>, true>();
 	} else if constexpr (particleType == ParticleType::StochasticStellarPop) {
 		return expandEnumNames<StochasticStellarPopParticleRealIdx, StochasticStellarPopParticleRealComps<problem_t>, true>();
+	} else if constexpr (particleType == ParticleType::Star) {
+		return expandEnumNames<StarParticleDataIdx, StarParticleRealComps<problem_t>, true>();
 	} else if constexpr (particleType == ParticleType::Sink) {
 		return expandEnumNames<SinkParticleRealIdx, SinkParticleRealComps, false>();
 	} else if constexpr (particleType == ParticleType::Test) {
@@ -407,6 +432,8 @@ template <ParticleType particleType, typename problem_t> auto getParticleIntComp
 								     // No integer components
 	} else if constexpr (particleType == ParticleType::CICRad) { // NOLINT
 								     // No integer components
+	} else if constexpr (particleType == ParticleType::Star) {   // NOLINT
+								     // No integer components
 	} else if constexpr (particleType == ParticleType::StochasticStellarPop) {
 		const std::vector<std::string> enum_names = amrex::getEnumNameStrings<StochasticStellarPopParticleIntIdx>();
 		names = {enum_names.begin(), enum_names.end()};
@@ -453,6 +480,15 @@ inline auto get_units_data() -> const auto &
 	       {"death_density", {1, -3, 0, 0}},
 	       {"mass_at_birth", {1, 0, 0, 0}},
 	       {"luminosity", {-1, 2, -3, 0}}}}},
+	    {ParticleType::Star,
+	     {{{"mass", {1, 0, 0, 0}},
+	       {"vx", {0, 1, -1, 0}},
+	       {"vy", {0, 1, -1, 0}},
+	       {"vz", {0, 1, -1, 0}},
+	       {"birth_time", {0, 0, 1, 0}},
+	       {"mdot", {1, 0, -1, 0}},
+	       {"radius", {0, 1, 0, 0}},
+	       {"luminosity", {-1, 2, -3, 0}}}}},
 	    {ParticleType::Sink,
 	     {{{"mass", {1, 0, 0, 0}},
 	       {"vx", {0, 1, -1, 0}},

src/particles/particle_update.hpp

@@ -6,6 +6,7 @@
 #include "particle_radiation.hpp"
 #include "particle_types.hpp"
 #include "physics_info.hpp"
+#include "starparticle_radiation.hpp"
 
 namespace quokka
 {
@@ -105,6 +106,30 @@ template <> struct ParticlePropertyUpdateTraits<ParticleType::StochasticStellarP
 	}
 };
 
+#if AMREX_SPACEDIM == 3
+// Specialization for Star particles: dispatches to the modular stellar-evolution framework.
+template <> struct ParticlePropertyUpdateTraits<ParticleType::Star> : ParticlePropertyUpdateBase<ParticleType::Star> {
+	template <typename problem_t, typename ParticleType, int Nout>
+	AMREX_GPU_DEVICE AMREX_FORCE_INLINE static void updateProperties(ParticleType &p, amrex::Real current_time, amrex::Real dt,
+									 LuminosityGpuConstTables<Nout> const &gpu_tables) noexcept
+	{
+		StellarUpdate::updateStellarProperties<problem_t>(p, current_time, dt, gpu_tables);
+	}
+
+	template <typename problem_t, typename ContainerType>
+	static void updateParticleProperties(ContainerType *container, amrex::Real current_time, amrex::Real dt) noexcept
+	{
+		const BL_PROFILE("ParticlePropertyUpdateTraits<Star>::updateParticleProperties()");
+		if (container == nullptr) {
+			return;
+		}
+		constexpr int nGroups = Physics_Traits<problem_t>::nGroups;
+		LuminosityGpuConstTables<nGroups> const gpu_tables{}; // unused by the stellar model, passed for signature parity
+		applyUpdate<problem_t, ContainerType>(container, current_time, dt, gpu_tables);
+	}
+};
+#endif // AMREX_SPACEDIM == 3
+
 // // Specialization for StochasticStellarPop particles from a simple analytical formula
 // // This is kept for debugging purpose.
 // template <> struct ParticlePropertyUpdateTraits<ParticleType::StochasticStellarPop> {

src/particles/star_particle_indices.H

@@ -0,0 +1,39 @@
+#ifndef STAR_PARTICLE_INDICES_H_
+#define STAR_PARTICLE_INDICES_H_
+
+#include "AMReX_Enum.H"
+
+namespace quokka
+{
+
+// Star particle real-component layout.
+// Every Star particle has StarParticleFixedComps fixed scalars followed by
+// nGroups luminosity slots and model-defined extras (nExtraReal / nExtraInt).
+AMREX_ENUM(StarParticleDataIdx, // NOLINT
+	   mass,		// Mass of the particle
+	   vx,			// Velocity x
+	   vy,			// Velocity y
+	   vz,			// Velocity z
+	   birth_time,		// Simulation time when the particle was created
+	   mdot,		// Current mass accretion rate (set by the accretion module)
+	   radius,		// Stellar radius (set by the stellar-evolution model)
+	   lum			// Base index for luminosity components (MUST be last; expands to lum_0, lum_1, ... for nGroups)
+);
+
+constexpr int StarParticleMassIdx = static_cast<int>(StarParticleDataIdx::mass);
+constexpr int StarParticleVxIdx = static_cast<int>(StarParticleDataIdx::vx);
+constexpr int StarParticleVyIdx = static_cast<int>(StarParticleDataIdx::vy);
+constexpr int StarParticleVzIdx = static_cast<int>(StarParticleDataIdx::vz);
+constexpr int StarParticleBirthTimeIdx = static_cast<int>(StarParticleDataIdx::birth_time);
+constexpr int StarParticleMdotIdx = static_cast<int>(StarParticleDataIdx::mdot);
+constexpr int StarParticleRadiusIdx = static_cast<int>(StarParticleDataIdx::radius);
+constexpr int StarParticleLumIdx = static_cast<int>(StarParticleDataIdx::lum);
+
+static_assert(static_cast<int>(StarParticleDataIdx::lum) == 7);
+
+// Number of fixed components before the variable-length sections (nGroups + extras).
+constexpr int StarParticleFixedComps = 7;
+
+} // namespace quokka
+
+#endif // STAR_PARTICLE_INDICES_H_

src/particles/starparticle_radiation.hpp

@@ -0,0 +1,35 @@
+#ifndef STARPARTICLE_RADIATION_HPP_
+#define STARPARTICLE_RADIATION_HPP_
+
+#include "AMReX_Extension.H"
+#include "AMReX_GpuQualifiers.H"
+#include "AMReX_REAL.H"
+
+#include "particles/particle_radiation.hpp" // LuminosityGpuConstTables
+#include "particles/particle_types.hpp"	    // StarParticle*Idx, ParticleType, Particle_Traits
+
+#if AMREX_SPACEDIM == 3
+
+namespace quokka
+{
+
+// Framework dispatcher for per-particle stellar-evolution updates.
+// Passes the full particle real-data array to the model so it can read and
+// modify any component (mass, radius, luminosity groups, etc.).
+class StellarUpdate
+{
+      public:
+	template <typename problem_t, typename ParticleType, int Nout>
+	AMREX_GPU_DEVICE AMREX_FORCE_INLINE static void updateStellarProperties(ParticleType &p, amrex::Real /*current_time*/, amrex::Real dt,
+										LuminosityGpuConstTables<Nout> const & /*gpu_tables*/) noexcept
+	{
+		using Model = typename Particle_Traits<problem_t>::stellar_model;
+		Model::evolve(&p.rdata(0), Nout, dt);
+	}
+};
+
+} // namespace quokka
+
+#endif // AMREX_SPACEDIM == 3
+
+#endif // STARPARTICLE_RADIATION_HPP_