diff --git a/src/particles/PhysicsParticles.hpp b/src/particles/PhysicsParticles.hpp index 318073af90..b4f7e97ec0 100644 --- a/src/particles/PhysicsParticles.hpp +++ b/src/particles/PhysicsParticles.hpp @@ -9,12 +9,12 @@ #include #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" @@ -115,7 +115,8 @@ class PhysicsParticleDescriptorBase [[nodiscard]] virtual auto computeMaxParticleSpeed(int lev) const -> amrex::ValLocPair = 0; // Methods that are implemented for some but not all particle types, so they cannot be pure virtual - virtual void depositSN(const amrex::Vector &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 }; @@ -568,7 +569,7 @@ template } // 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()."); @@ -636,19 +637,21 @@ class StarParticleDescriptor : public PhysicsParticleDescriptor &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::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::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(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 @@ -678,6 +681,17 @@ template 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 { @@ -782,13 +796,12 @@ template class PhysicsParticleRegister } // Deposit supernova energy and momentum from all particles - void depositSN(const amrex::Vector &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); } } } diff --git a/src/particles/particle_deposition.hpp b/src/particles/particle_deposition.hpp index 42d59cbf37..0d9430723d 100644 --- a/src/particles/particle_deposition.hpp +++ b/src/particles/particle_deposition.hpp @@ -5,13 +5,16 @@ #include "AMReX_Extension.H" #include "AMReX_MultiFab.H" #include "AMReX_ParticleInterpolators.H" -#include "AMReX_ParticleMesh.H" #include "AMReX_REAL.H" +#include "hydro/hydro_system.hpp" #include "particles/particle_types.hpp" namespace quokka { +constexpr double cloudy_H_mass_fraction = 1.0 / (1.0 + 0.1 * 3.971); +constexpr double m_u = C::m_u; + //-------------------- Radiation depositions -------------------- // Functor for depositing radiation energy from particles onto the grid @@ -66,82 +69,355 @@ struct MassDeposition { }; //-------------------- Supernova depositions -------------------- +template +void SNDeposition(ContainerType *container, amrex::MultiFab &state, amrex::MultiFab &state_buffer, int lev, amrex::Real time, amrex::Real dt, int mass_index, + int evolutionStageIndex, int birthTimeIndex) +{ + constexpr int stencil_size = 3; + constexpr amrex::Real stencil_volume = 4.0 / 3.0 * M_PI * stencil_size * stencil_size * stencil_size; + constexpr amrex::Real stencil_weights[4][4][4] = // NOLINT + {{{0.00884198143074, 0.00884198143074, 0.00884198143074, 0.00416240696843}, + {0.00884198143074, 0.00884198143074, 0.00884198143074, 0.00262865918549}, + {0.00884198143074, 0.00884198143074, 0.00596795726055, 0.00005052308190}, + {0.00416240696843, 0.00262865918549, 0.00005052308190, 0.00000000000000}}, + {{0.00884198143074, 0.00884198143074, 0.00884198143074, 0.00262865918549}, + {0.00884198143074, 0.00884198143074, 0.00861063982859, 0.00119306623841}, + {0.00884198143074, 0.00861063982859, 0.00400459528385, 0.00000136166514}, + {0.00262865918549, 0.00119306623841, 0.00000136166514, 0.00000000000000}}, + {{0.00884198143074, 0.00884198143074, 0.00596795726055, 0.00005052308190}, + {0.00884198143074, 0.00861063982859, 0.00400459528385, 0.00000136166514}, + {0.00596795726055, 0.00400459528385, 0.00045652034325, 0.00000000000000}, + {0.00005052308190, 0.00000136166514, 0.00000000000000, 0.00000000000000}}, + {{0.00416240696843, 0.00262865918549, 0.00005052308190, 0.00000000000000}, + {0.00262865918549, 0.00119306623841, 0.00000136166514, 0.00000000000000}, + {0.00005052308190, 0.00000136166514, 0.00000000000000, 0.00000000000000}, + {0.00000000000000, 0.00000000000000, 0.00000000000000, 0.00000000000000}}}; + + static_assert(stencil_size <= 3, + "stencil_size must be <= 3"); // stencil_size must be <= n_ghost - 1 = 3. SN particle may drift 1 cell before being deposited. + + // copy host variables to device + const SNScheme SN_scheme_d = SN_scheme; + + const amrex::Real step_end_time = time + dt; + + constexpr double E_blast = 1.0e51; // ergs + constexpr double m_ej = 10.0 * C::M_solar; // ejecta mass in cgs + constexpr double m_dead_min = 1.4 * C::M_solar; // minimum mass of a dead star + constexpr double p_snr_0 = 2.8e5 * C::M_solar * 1.0e5; // SN terminal momentum in cgs + + // Zero the buffer for each particle type + state_buffer.setVal(0.0); + + // Step 1: Local deposition within each box + for (typename ContainerType::ParIterType pti(*container, lev); pti.isValid(); ++pti) { + // Get the particle array of structs + auto &particles = pti.GetArrayOfStructs(); + auto *pData = particles().data(); + const amrex::Long np = pti.numParticles(); + + // Get the local deposit array for this box + const auto &local_buffer = state_buffer.array(pti); + const auto &local_state = state.array(pti); + + // Get geometry information for this level + const auto &geom = container->Geom(lev); + const auto plo = geom.ProbLoArray(); + const auto dxi = geom.InvCellSizeArray(); + const auto dx = geom.CellSizeArray(); + + // Calculate inverse cell volume + const amrex::Real vol_inverse = AMREX_D_TERM(dxi[0], *dxi[1], *dxi[2]); + const amrex::Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); + + // Deposit particle data into the local buffer + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int64_t idx) { + auto &p = pData[idx]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic) -// Functor for depositing supernova energy and momentum from particles onto the grid -// This is a simplified version of the SNDeposition functor that deposits mass and energy uniformly -// to 5³ cells centered on the particle's cell. It is used for testing purposes. -// Note: the deposition radius must be <= nghost_cc_ -struct SNDeposition { - double step_end_time{}; // Current simulation time - int start_part_comp{}; // Starting component in particle data - int start_mesh_comp{}; // Starting component in mesh data - int birthTimeIndex{}; // Index for particle birth time - int evolutionStageIndex{}; // Index for particle evolution stage - double SN_time = particle_param2; - - // For some unknown reason, stencil_width < 3 results in larger error in SNR mass when a particle is at the refinement boundary. - static constexpr int stencil_width = 3; - - // Abort if stencil_width > nghost_cc_ - 1. - // The particle can drift one cell out of the valid zone during kickParticlesAllLevels(). - // A stencil_width > nghost_cc_ - 1 would result in particles depositing energy/momentum outside the ghost zones. - // We can't use AMRSimulation::nghost_cc_ and have to hard-code 3 here because we don't have a problem_t template parameter. - static_assert(stencil_width <= 3, "stencil_width must be <= nghost_cc_"); - - // Operator to perform supernova deposition using cloud-in-cell approach - template - AMREX_GPU_DEVICE AMREX_FORCE_INLINE void operator()(const ContainerType &p, amrex::Array4 const &state, - amrex::GpuArray const &plo, - amrex::GpuArray const &dxi) const noexcept - { - // Check if the particle has an integer component for evolution stage - if constexpr (ContainerType::NInt > 0) { // Check if this is a supernova progenitor - bool is_sn_progenitor = false; - if (evolutionStageIndex >= 0) { - is_sn_progenitor = (p.idata(evolutionStageIndex) == static_cast(StellarEvolutionStage::SNProgenitor)); - } + const bool is_sn_progenitor = (p.idata(evolutionStageIndex) == static_cast(StellarEvolutionStage::SNProgenitor)); + + if (is_sn_progenitor && step_end_time > p.rdata(birthTimeIndex + 1)) { + // Update the particle's evolution stage to SNRemnant + p.idata(evolutionStageIndex) = static_cast(StellarEvolutionStage::SNRemnant); + + // update the particle's mass: subtract ejecta mass + const amrex::Real mass_dead_star = p.rdata(mass_index) - m_ej; + // AMREX_ASSERT_WITH_MESSAGE(mass_dead_star > 0.0, "SN progenitor mass should be greater than ejecta mass (10 M_sun)"); + p.rdata(mass_index) = std::max(m_dead_min, mass_dead_star); + + // get particle velocity and kinetic energy + const amrex::Real p_vx = p.rdata(mass_index + 1); + const amrex::Real p_vy = p.rdata(mass_index + 2); + const amrex::Real p_vz = p.rdata(mass_index + 3); + const amrex::Real SN_kin_energy = 0.5 * m_ej * (p_vx * p_vx + p_vy * p_vy + p_vz * p_vz); - if (is_sn_progenitor && step_end_time > p.rdata(birthTimeIndex) + SN_time) { // Find the cell containing the particle - int base_i = static_cast(amrex::Math::floor((p.pos(0) - plo[0]) * dxi[0])); - int base_j = static_cast(amrex::Math::floor((p.pos(1) - plo[1]) * dxi[1])); - int base_k = static_cast(amrex::Math::floor((p.pos(2) - plo[2]) * dxi[2])); - - // Calculate the volume factor for normalization (5³ cells) - const int num_cells = (2 * stencil_width + 1) * (2 * stencil_width + 1) * (2 * stencil_width + 1); - const amrex::Real vol_factor = (AMREX_D_TERM(dxi[0], *dxi[1], *dxi[2])) / num_cells; - - // Deposit evenly to 5³ cells centered on the particle's cell - const amrex::Real pdensity = p.rdata(start_part_comp) * vol_factor; - const amrex::Real penergy = pdensity; // for testing: energy = mass - const amrex::Real pmomentum = 0.0; // for testing: momentum = 0 - - for (int kk = -stencil_width; kk <= stencil_width; ++kk) { - for (int jj = -stencil_width; jj <= stencil_width; ++jj) { - for (int ii = -stencil_width; ii <= stencil_width; ++ii) { - // Add the contribution to each cell - // We assume start_mesh_comp is the density index followed by the momentum indices and then the energy - // index - amrex::Gpu::Atomic::AddNoRet(&state(base_i + ii, base_j + jj, base_k + kk, start_mesh_comp), pdensity); - amrex::Gpu::Atomic::AddNoRet(&state(base_i + ii, base_j + jj, base_k + kk, start_mesh_comp + 1), - pmomentum); - amrex::Gpu::Atomic::AddNoRet(&state(base_i + ii, base_j + jj, base_k + kk, start_mesh_comp + 2), - pmomentum); - amrex::Gpu::Atomic::AddNoRet(&state(base_i + ii, base_j + jj, base_k + kk, start_mesh_comp + 3), - pmomentum); - amrex::Gpu::Atomic::AddNoRet(&state(base_i + ii, base_j + jj, base_k + kk, start_mesh_comp + 4), - penergy); + int ix = static_cast(amrex::Math::floor((p.pos(0) - plo[0]) * dxi[0])); + int iy = static_cast(amrex::Math::floor((p.pos(1) - plo[1]) * dxi[1])); + int iz = static_cast(amrex::Math::floor((p.pos(2) - plo[2]) * dxi[2])); + + amrex::Real avg_density = 0.0; + for (int ii = -stencil_size; ii <= stencil_size; ++ii) { + for (int jj = -stencil_size; jj <= stencil_size; ++jj) { + for (int kk = -stencil_size; kk <= stencil_size; ++kk) { + const int iii = std::abs(ii); + const int jjj = std::abs(jj); + const int kkk = std::abs(kk); + const double kernel = stencil_weights[iii][jjj][kkk]; + avg_density += kernel * local_state(ix + ii, iy + jj, iz + kk, HydroSystem::density_index); } } } - // Note: We cannot modify the particle here because it's passed as const reference - // The evolution stage update is now handled by the updateEvolutionStage function + if (SN_scheme_d == SNScheme::SN_thermal_only) { + // Deposit mass and energy into (2 * stencil_width + 1)³ cells centered on the particle's cell + for (int ii = -stencil_size; ii <= stencil_size; ++ii) { + for (int jj = -stencil_size; jj <= stencil_size; ++jj) { + for (int kk = -stencil_size; kk <= stencil_size; ++kk) { + const int iii = std::abs(ii); + const int jjj = std::abs(jj); + const int kkk = std::abs(kk); + const double kernel_times_vol_inverse = stencil_weights[iii][jjj][kkk] * vol_inverse; + const amrex::Real SNR_rho_per_cell = m_ej * kernel_times_vol_inverse; + const amrex::Real SNR_energy_per_cell = (E_blast + SN_kin_energy) * kernel_times_vol_inverse; + const amrex::Real SNR_px_per_cell = m_ej * p_vx * kernel_times_vol_inverse; + const amrex::Real SNR_py_per_cell = m_ej * p_vy * kernel_times_vol_inverse; + const amrex::Real SNR_pz_per_cell = m_ej * p_vz * kernel_times_vol_inverse; + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ix + ii, iy + jj, iz + kk, HydroSystem::density_index), + SNR_rho_per_cell); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ix + ii, iy + jj, iz + kk, HydroSystem::x1Momentum_index), + SNR_px_per_cell); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ix + ii, iy + jj, iz + kk, HydroSystem::x2Momentum_index), + SNR_py_per_cell); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ix + ii, iy + jj, iz + kk, HydroSystem::x3Momentum_index), + SNR_pz_per_cell); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ix + ii, iy + jj, iz + kk, HydroSystem::energy_index), + SNR_energy_per_cell); + } + } + } + } else { + const double n_H_amb = avg_density * cloudy_H_mass_fraction / m_u; + const amrex::Real M_snr = (avg_density * stencil_volume * vol) + m_ej; // SNR mass + const amrex::Real M_sf = 1679.0 * C::M_solar * std::pow(n_H_amb, -0.26); // Shell-formation mass + const amrex::Real RM = M_snr / M_sf; // R_M factor = M_snr / M_sf + const amrex::Real p_snr = p_snr_0 * std::pow(n_H_amb, -0.17); // = 1.89e5 when n = 10 + + // fraction of terminal SN momentum to go to gas momentum + amrex::Real f_factor = 1.0; + if (SN_scheme_d == SNScheme::SN_thermal_or_thermal_momentum) { + if (RM < 1.0) { + if (RM > 0.027) { + f_factor = 0.529 * std::sqrt(RM); // f^2 = 0.28. 28% kinetic (Kim & Ostriker 2017) + } else { + f_factor = 0.0; // pure thermal in well-resolved limit + } + } + } else if (SN_scheme_d == SNScheme::SN_thermal_kinetic_or_thermal_momentum) { + if (RM < 1.0) { + if (RM > 0.027) { + f_factor = 0.529 * std::sqrt(RM); // f^2 = 0.28. 28% kinetic (Kim & Ostriker 2017) + } else { + f_factor = 1.414 * std::sqrt(RM); // pure kinetic in well-resolved limit: (f^2 / RM) * (p_snr^2 + // / 2 M_sf) = 2 * (p_snr^2 / 2 M_sf) ~= 1e51 erg + } + } + } + // SNScheme::SN_pure_kinetic_or_thermal_momentum: keep f_factor = 1.0 + + // // log RM and f_factor, for debugging on CPU. + // if (particle_verbose) { + // printf("SNR logging -- RM: %.2e, f_factor: %.2e\n", RM, f_factor); + // } + + for (int ii = ix - stencil_size; ii <= ix + stencil_size; ++ii) { + for (int jj = iy - stencil_size; jj <= iy + stencil_size; ++jj) { + for (int kk = iz - stencil_size; kk <= iz + stencil_size; ++kk) { + const int iii = std::abs(ii - ix); + const int jjj = std::abs(jj - iy); + const int kkk = std::abs(kk - iz); + const double kernel_times_vol_inverse = stencil_weights[iii][jjj][kkk] * vol_inverse; + + const double delta_x = (ii + 0.5) * dx[0] + plo[0] - p.pos(0); + const double delta_y = (jj + 0.5) * dx[1] + plo[1] - p.pos(1); + const double delta_z = (kk + 0.5) * dx[2] + plo[2] - p.pos(2); + const double r_sq = (delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z); + + const amrex::Real delta_rho_i = m_ej * kernel_times_vol_inverse; + const amrex::Real e_snr_per_cell = (E_blast + SN_kin_energy) * kernel_times_vol_inverse; + const amrex::Real momentum_per_cell = f_factor * p_snr * kernel_times_vol_inverse; + + // Compute unit vector from particle to cell center + const amrex::Real r = std::sqrt(r_sq); + // Avoid division by zero + const amrex::Real inv_r = (r > 0.0) ? 1.0 / r : 0.0; + + // unit vector from particle to cell center + const amrex::Real r_hat_x = delta_x * inv_r; + const amrex::Real r_hat_y = delta_y * inv_r; + const amrex::Real r_hat_z = delta_z * inv_r; + + const double rho = local_state(ii, jj, kk, HydroSystem::density_index); + const double px = local_state(ii, jj, kk, HydroSystem::x1Momentum_index); + const double py = local_state(ii, jj, kk, HydroSystem::x2Momentum_index); + const double pz = local_state(ii, jj, kk, HydroSystem::x3Momentum_index); + + // Compute momentum directed along unit vector + const double dpx = (delta_rho_i * px / rho) + (momentum_per_cell * r_hat_x); + const double dpy = (delta_rho_i * py / rho) + (momentum_per_cell * r_hat_y); + const double dpz = (delta_rho_i * pz / rho) + (momentum_per_cell * r_hat_z); + + amrex::Gpu::Atomic::AddNoRet(&local_buffer(ii, jj, kk, HydroSystem::density_index), + delta_rho_i); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ii, jj, kk, HydroSystem::x1Momentum_index), dpx); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ii, jj, kk, HydroSystem::x2Momentum_index), dpy); + amrex::Gpu::Atomic::AddNoRet( + &local_buffer(ii, jj, kk, HydroSystem::x3Momentum_index), dpz); + amrex::Gpu::Atomic::AddNoRet(&local_buffer(ii, jj, kk, HydroSystem::energy_index), + e_snr_per_cell); + } + } + } + } } - } + }); } -}; + + // Step 2: Sum boundary values + state_buffer.SumBoundary(container->Geom(lev).periodicity()); + + // Step 3: Add the buffer to the state + for (amrex::MFIter mfi(state); mfi.isValid(); ++mfi) { + const amrex::Box &box = mfi.validbox(); + auto const &local_state = state.array(mfi); + auto const &local_buffer = state_buffer.array(mfi); + + // add buffer to state + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + if (SN_scheme_d == SNScheme::SN_thermal_only) { + // For SN_thermal_only, the buffer contains only mass and energy (and a small amount of momentum), so it's safe to add the + // buffer directly to the state. + const double rho_new = + local_state(i, j, k, HydroSystem::density_index) + local_buffer(i, j, k, HydroSystem::density_index); + const double px_new = local_state(i, j, k, HydroSystem::x1Momentum_index) + + local_buffer(i, j, k, HydroSystem::x1Momentum_index); + const double py_new = local_state(i, j, k, HydroSystem::x2Momentum_index) + + local_buffer(i, j, k, HydroSystem::x2Momentum_index); + const double pz_new = local_state(i, j, k, HydroSystem::x3Momentum_index) + + local_buffer(i, j, k, HydroSystem::x3Momentum_index); + const double e_new = + local_state(i, j, k, HydroSystem::energy_index) + local_buffer(i, j, k, HydroSystem::energy_index); + const double e_int_new = e_new - (0.5 * ((px_new * px_new) + (py_new * py_new) + (pz_new * pz_new)) / rho_new); + + local_state(i, j, k, HydroSystem::density_index) = rho_new; + local_state(i, j, k, HydroSystem::x1Momentum_index) = px_new; + local_state(i, j, k, HydroSystem::x2Momentum_index) = py_new; + local_state(i, j, k, HydroSystem::x3Momentum_index) = pz_new; + local_state(i, j, k, HydroSystem::internalEnergy_index) = e_int_new; + local_state(i, j, k, HydroSystem::energy_index) = e_new; + } else { + // For SN_thermal_or_thermal_momentum, SN_thermal_kinetic_or_thermal_momentum, and SN_pure_kinetic_or_thermal_momentum, + // the buffer contains mass, momentum, and energy. We need to add the buffer to the state in a way that guarantees + // that the internal energy is positive. In fact, we demand that the cell temperature should not decrease. + const double rho = local_state(i, j, k, HydroSystem::density_index); + const double px = local_state(i, j, k, HydroSystem::x1Momentum_index); + const double py = local_state(i, j, k, HydroSystem::x2Momentum_index); + const double pz = local_state(i, j, k, HydroSystem::x3Momentum_index); + const double e_int = local_state(i, j, k, HydroSystem::internalEnergy_index); + const double e_tot = local_state(i, j, k, HydroSystem::energy_index); + + const double d_rho = local_buffer(i, j, k, HydroSystem::density_index); + const double d_px = local_buffer(i, j, k, HydroSystem::x1Momentum_index); + const double d_py = local_buffer(i, j, k, HydroSystem::x2Momentum_index); + const double d_pz = local_buffer(i, j, k, HydroSystem::x3Momentum_index); + const double d_e = local_buffer(i, j, k, HydroSystem::energy_index); + + const double rho_new = rho + d_rho; + double px_new = px + d_px; + double py_new = py + d_py; + double pz_new = pz + d_pz; + double e_tot_new = e_tot + d_e; + + const double d_e_int_d_rho = e_int / rho; + const double e_int_new_tmp = d_e_int_d_rho * rho_new; + const double e_int_plus_kinetic = e_int_new_tmp + (0.5 * ((px_new * px_new) + (py_new * py_new) + (pz_new * pz_new)) / rho_new); + + if (e_int_plus_kinetic <= e_tot_new * (1.0 + 1.0e-10)) { + e_tot_new = std::max(e_int_plus_kinetic, e_tot_new); + } else { + // find the lambda such that e_int_plus_kinetic == e_tot_new + const double e_kinetic_max = e_tot_new - e_int_new_tmp; + AMREX_ASSERT(e_kinetic_max >= 0.0); + + // If this assertion fails, it means the SN energy (10^51 erg) is not enough to accelerate the + // SNR to match the velocity of the gas. This may ONLY happen when the SN star is nearly static + // at first and the background gas is moving at a speed >3171 km/s. This should NOT happen in + // practice. + // TODO(cch): deal with the special case where this assertion fails. + AMREX_ASSERT(0.5 * (px * px + py * py + pz * pz) / rho_new <= e_kinetic_max); + + // Find analytical solution of the following equation: + // d_p^2 lambda^2 + 2 d_p p lambda + p^2 - 2 rho_new e_kinetic_max = 0 + // This quadratic equation, denoted as F(x) = 0, has some known properties which make the + // solution well constrained: + // 1. F(0) < 0 (see assertion above) + // 2. F(1) > 0 (otherwise we won't be in this else clause) + // Therefore, there must be one and only one solution in the range [0, 1], and it's the bigger + // one of the two solutions (so we take the plus sign in the quadratic formula). + const double a = (d_px * d_px) + (d_py * d_py) + (d_pz * d_pz); + const double b = 2.0 * (px * d_px + py * d_py + pz * d_pz); + const double c = (px * px) + (py * py) + (pz * pz) - (2.0 * rho_new * e_kinetic_max); + const double discriminant = (b * b) - (4.0 * a * c); + AMREX_ASSERT(discriminant >= 0.0); + const double lambda = (-b + std::sqrt(discriminant)) / (2.0 * a); + AMREX_ASSERT(lambda >= 0.0); + AMREX_ASSERT(lambda <= 1.0); + + // assert that lambda is a valid solution + AMREX_ASSERT_WITH_MESSAGE( + std::abs((0.5 * + ((px + lambda * d_px) * (px + lambda * d_px) + (py + lambda * d_py) * (py + lambda * d_py) + + (pz + lambda * d_pz) * (pz + lambda * d_pz)) / + rho_new) - + e_kinetic_max) <= 1.0e-10 * e_kinetic_max, + "lambda is not a valid solution. This should NOT happen. @chongchonghe should be responsible for this."); + + px_new = px + lambda * d_px; + py_new = py + lambda * d_py; + pz_new = pz + lambda * d_pz; + } + + const double e_int_new = e_tot_new - (0.5 * ((px_new * px_new) + (py_new * py_new) + (pz_new * pz_new)) / rho_new); + local_state(i, j, k, HydroSystem::density_index) = rho_new; + local_state(i, j, k, HydroSystem::x1Momentum_index) = px_new; + local_state(i, j, k, HydroSystem::x2Momentum_index) = py_new; + local_state(i, j, k, HydroSystem::x3Momentum_index) = pz_new; + local_state(i, j, k, HydroSystem::internalEnergy_index) = e_int_new; + local_state(i, j, k, HydroSystem::energy_index) = e_tot_new; + + // // log the state, for debugging on CPU. + // if (d_rho / rho > 1.0e-12) { + // // print original rho, px, py, pz, e_int, e_tot; new rho_new, px_new, py_new, pz_new, e_int_new, + // // e_tot_new + // printf("original: rho = %e, px = %e, py = %e, pz = %e, e_int = %e, e_tot = %e\n", rho, px, py, pz, e_int, e_tot); + // printf("new: rho_new = %e, px_new = %e, py_new = %e, pz_new = %e, e_int_new = %e, e_tot_new = %e\n", rho_new, px_new, + // py_new, pz_new, e_int_new, e_tot_new); + // // print d e_int / d e_tot + // printf("d e_int / d e_tot = %e\n", (e_int_new - e_int) / (e_tot_new - e_tot)); + // printf("e_int / rho = %e, e_int_new / rho_new = %e\n", e_int / rho, e_int_new / rho_new); + // } + } + }); + } +} // Function to update particle evolution stages from SNProgenitor to SNRemnant template @@ -163,8 +439,7 @@ void updateEvolutionStage(ContainerType *container, int lev_min, amrex::Real ste // Check if this is a supernova progenitor const bool is_sn_progenitor = (p.idata(evolutionStageIndex) == static_cast(StellarEvolutionStage::SNProgenitor)); - // Update the particle's evolution stage to SNRemnant if it is a progenitor and the end of the current time step is greater than - // the death time. + // Update the particle's evolution stage to SNRemnant if it's time if (is_sn_progenitor && step_end_time > p.rdata(birthTimeIndex + 1)) { p.idata(evolutionStageIndex) = static_cast(StellarEvolutionStage::SNRemnant); } diff --git a/src/particles/particle_types.hpp b/src/particles/particle_types.hpp index 17802dc9f9..e818b60ab3 100644 --- a/src/particles/particle_types.hpp +++ b/src/particles/particle_types.hpp @@ -75,15 +75,13 @@ enum class ParticleType { Test // Test particles with all features enabled }; -// Global particle parameters -// The 'inline' keyword is used here to avoid multiple definition errors when this header -// is included in multiple source files. It ensures that all translation units that include -// this header will refer to the same instance of these variables, rather than creating -// 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 +// Enum for SN schemes: ThermalOnly, ThermalAndMomentum +enum class SNScheme { + SN_thermal_only, // pure thermal + SN_thermal_or_thermal_momentum, // pure thermal (RM<1) or thermal+momentum (RM>=1) + SN_thermal_kinetic_or_thermal_momentum, // thermal+kinetic (RM<1) or thermal+momentum (RM>=1) + SN_pure_kinetic_or_thermal_momentum // pure kinetic (RM<1) or thermal+momentum (RM>=1) +}; //-------------------- Radiation particles -------------------- @@ -270,6 +268,26 @@ inline auto get_units_data() -> const auto & // 1. For massive particles, velocity components start after mass // 2. Birth time, if existing, is always followed by death time +// Global particle parameters +// The 'inline' keyword is used here to avoid multiple definition errors when this header +// is included in multiple source files. It ensures that all translation units that include +// this header will refer to the same instance of these variables, rather than creating +// their own copies. + +// Disable SN feedback when a particle evolves from SNProgenitor to SNRemnant +inline bool disable_SN_feedback = false; // NOLINT + +// Placeholder parameters for particles. Used in gravity_3d.cpp tests +inline amrex::Real particle_param1 = -1.0; // NOLINT +inline amrex::Real particle_param2 = -1.0; // NOLINT +inline amrex::Real particle_param3 = -1.0; // NOLINT + +// Scheme for SN feedback +inline SNScheme SN_scheme = SNScheme::SN_thermal_only; // NOLINT + +// Verbosity for particle operations +inline int particle_verbose = 0; // NOLINT print particle logistics + // Function to parse particle parameters from input file // The 'inline' keyword allows this function to be defined in a header file without // causing multiple definition errors when the header is included in multiple source files. @@ -280,9 +298,17 @@ inline void particleParmParse() { // Parse particle parameters const amrex::ParmParse pp("particles"); + pp.query("disable_SN_feedback", disable_SN_feedback); pp.query("param1", particle_param1); pp.query("param2", particle_param2); pp.query("param3", particle_param3); + + // Handle SNScheme enum + int sn_scheme_int = static_cast(SN_scheme); + pp.query("SN_scheme", sn_scheme_int); + SN_scheme = static_cast(sn_scheme_int); + + // Handle integer verbose flag pp.query("verbose", particle_verbose); } diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index d447dd7f16..4e9838f0a0 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -3,6 +3,7 @@ add_subdirectory(Advection2D) add_subdirectory(AdvectionSemiellipse) add_subdirectory(ParticleCreation) +add_subdirectory(SN) add_subdirectory(HydroBlast2D) add_subdirectory(HydroBlast3D) @@ -68,4 +69,4 @@ add_subdirectory(PrimordialChem) add_subdirectory(PopIII) add_subdirectory(ShockCloud) add_subdirectory(StarCluster) -add_subdirectory(SphericalCollapse) \ No newline at end of file +add_subdirectory(SphericalCollapse) diff --git a/src/problems/SN/CMakeLists.txt b/src/problems/SN/CMakeLists.txt new file mode 100644 index 0000000000..4ff1babc22 --- /dev/null +++ b/src/problems/SN/CMakeLists.txt @@ -0,0 +1,8 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(test_SN test_SN.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_SN) + endif() + + add_test(NAME SN COMMAND test_SN SN.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/problems/SN/test_SN.cpp b/src/problems/SN/test_SN.cpp new file mode 100644 index 0000000000..ee88294c85 --- /dev/null +++ b/src/problems/SN/test_SN.cpp @@ -0,0 +1,226 @@ +/// \file test_SN.cpp +/// \brief Defines a test problem for supernova feedback. +/// + +#include "AMReX.H" +#include "AMReX_BC_TYPES.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParmParse.H" +#include "AMReX_Print.H" +#include "AMReX_SPACE.H" + +#include "QuokkaSimulation.hpp" +#include "fundamental_constants.H" +#include "hydro/hydro_system.hpp" +#include "particles/particle_types.hpp" + +struct SNProblem { +}; + +static bool refine_half_domain = false; // NOLINT + +static double max_Eint_global = 0.0; // NOLINT + +static std::string SN_particles_file = "SN_particles.txt"; // NOLINT + +constexpr double mu = 1.0 * C::m_u; +// constexpr double mu = 1.295 * C::m_u; // neutral gas +constexpr double gamma_ = 5. / 3.; +const double CV = 1. / (gamma_ - 1.) / mu * C::k_B; +const double cloudy_H_mass_fraction = 1.0 / (1.0 + 0.1 * 3.971); +const double year = 3.15576e+07; // in seconds +const double mass_SNR = 10.0 * C::M_solar; +const int n_SNR = 2; + +static double n_amb = 1.0; // ambient density (g cm^-3) // NOLINT +static double T_amb = 100.0; // ambient temperature (K) // NOLINT +static double t_stop = 3.0e5; // stop time (yr) // NOLINT + +template <> struct Particle_Traits { + // static constexpr ParticleSwitch particle_switch = ParticleSwitch::None; + static constexpr ParticleSwitch particle_switch = ParticleSwitch::StochasticStellarPop; +}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = gamma_; + static constexpr double mean_molecular_weight = mu; +}; + +template <> struct HydroSystem_Traits { + static constexpr bool reconstruct_eint = true; // need to reconstruct temperature +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = false; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = 1; // number of radiation groups + static constexpr UnitSystem unit_system = UnitSystem::CGS; +}; + +template <> void QuokkaSimulation::createInitialStochasticStellarPopParticles() +{ + // read particles from ASCII file + const int nreal_extra = 7; // mass vx vy vz birth_time death_time lum + StochasticStellarPopParticles->SetVerbose(1); + StochasticStellarPopParticles->InitFromAsciiFile(SN_particles_file, nreal_extra, nullptr); + + // Loop over all particle at all levels and set first integer component to SNProgenitor + for (int lev = 0; lev <= StochasticStellarPopParticles->finestLevel(); ++lev) { + auto &particles = StochasticStellarPopParticles->GetParticles(lev); + + for (auto &kv : particles) { + auto &particle_array = kv.second.GetArrayOfStructs(); + const int np = particle_array.numParticles(); + auto *pdata = particle_array().data(); + + // Launch GPU kernel to set integer components + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int i) { + auto &p = pdata[i]; // NOLINT + p.idata(0) = static_cast(quokka::StellarEvolutionStage::SNProgenitor); + }); + } + } + + // Ensure GPU operations are complete + amrex::Gpu::streamSynchronize(); +} + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + const double rho_bg = n_amb * C::m_u / cloudy_H_mass_fraction; + const double E0 = CV * T_amb * rho_bg; + const double rho = rho_bg; + const double rho_e = E0; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + state_cc(i, j, k, HydroSystem::density_index) = rho; + state_cc(i, j, k, HydroSystem::x1Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::energy_index) = rho_e; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = rho_e; + }); +} + +template <> void QuokkaSimulation::ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real /*time*/, int /*ngrow*/) +{ + // tag cells for refinement: static mesh refinement for the whole domain (if refine_half_domain is false) or for x > 0 (if refine_half_domain is true) + + auto const &dx = geom[lev].CellSizeArray(); + auto const &plo = geom[lev].ProbLoArray(); + auto const &phi = geom[lev].ProbHiArray(); + const bool refine_half_domain_ = refine_half_domain; + + for (amrex::MFIter mfi(state_new_cc_[lev]); mfi.isValid(); ++mfi) { + const amrex::Box &box = mfi.validbox(); + const auto tag = tags.array(mfi); + + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const double x_frac = ((i + 0.5) * dx[0]) / (phi[0] - plo[0]); + const double y_frac = ((j + 0.5) * dx[1]) / (phi[1] - plo[1]); + const double z_frac = ((k + 0.5) * dx[2]) / (phi[2] - plo[2]); + if (!refine_half_domain_ || (x_frac >= 0.7 && x_frac <= 0.8 && y_frac >= 0.3 && y_frac <= 0.7 && z_frac >= 0.3 && z_frac <= 0.7)) { + tag(i, j, k) = amrex::TagBox::SET; + } + }); + } +} + +template <> void QuokkaSimulation::computeAfterTimestep() +{ + // find the maximum temperature in the state_new_cc_[0] + const double max_internal_energy_density = state_new_cc_[0].max(HydroSystem::internalEnergy_index); + max_Eint_global = std::max(max_Eint_global, max_internal_energy_density); +} + +auto problem_main() -> int +{ + auto isNormalComp = [=](int n, int dim) { + if ((n == HydroSystem::x1Momentum_index) && (dim == 0)) { + return true; + } + if ((n == HydroSystem::x2Momentum_index) && (dim == 1)) { + return true; + } + if ((n == HydroSystem::x3Momentum_index) && (dim == 2)) { + return true; + } + return false; + }; + + const int ncomp_cc = Physics_Indices::nvarTotal_cc; + amrex::Vector BCs_cc(ncomp_cc); + for (int n = 0; n < ncomp_cc; ++n) { + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + // // periodic boundaries + // BCs_cc[n].setLo(i, amrex::BCType::int_dir); + // BCs_cc[n].setHi(i, amrex::BCType::int_dir); + // octant symmetry + if (isNormalComp(n, i)) { + BCs_cc[n].setLo(i, amrex::BCType::reflect_odd); + BCs_cc[n].setHi(i, amrex::BCType::reflect_odd); + } else { + BCs_cc[n].setLo(i, amrex::BCType::reflect_even); + BCs_cc[n].setHi(i, amrex::BCType::reflect_even); + } + } + } + + // get n_amb from the input file + amrex::ParmParse const pp("problem"); + pp.query("n_amb", n_amb); + pp.query("T_amb", T_amb); + pp.query("t_stop", t_stop); + pp.query("SN_particles_file", SN_particles_file); + pp.query("refine_half_domain", refine_half_domain); + + // Problem initialization + QuokkaSimulation sim(BCs_cc); + + sim.reconstructionOrder_ = 3; // 2=PLM, 3=PPM + sim.stopTime_ = t_stop * year; + sim.cflNumber_ = 0.3; // *must* be less than 1/3 in 3D! + sim.initDt_ = 0.001 * year; + + // initialize + sim.setInitialConditions(); + + amrex::GpuArray const &dx0 = sim.geom[0].CellSizeArray(); + amrex::Real const vol = AMREX_D_TERM(dx0[0], *dx0[1], *dx0[2]); + amrex::Real const total_mass_init = sim.state_new_cc_[0].sum(HydroSystem::density_index) * vol; + + // evolve + sim.evolve(); + + amrex::Real const total_mass_final = sim.state_new_cc_[0].sum(HydroSystem::density_index) * vol; + const amrex::Real mass_increase = total_mass_final - total_mass_init; + amrex::Print() << "----------------- Problem diagnostics -----------------" << "\n"; + amrex::Print() << "Total mass increase: " << mass_increase << "\n"; + amrex::Print() << "Expected total mass increase: " << n_SNR * mass_SNR << "\n"; + const double mass_increase_rel_err = std::abs(mass_increase - n_SNR * mass_SNR) / (n_SNR * mass_SNR); + amrex::Print() << "Mass increase relative error: " << mass_increase_rel_err << "\n"; + const double mass_increase_rel_err_tol = 1.0e-10; + + const amrex::Real max_internal_energy = max_Eint_global * vol; + const amrex::Real expected_minimum_max_internal_energy = 1.0e51 / (7 * 7 * 7); // 1e51 erg energy into (2 * 3 + 1)^3 cells + int status = 1; + if (max_internal_energy > expected_minimum_max_internal_energy && mass_increase_rel_err < mass_increase_rel_err_tol) { + status = 0; + amrex::Print() << "Test passed. Max internal energy in cells: " << max_internal_energy << "\n"; + } else { + status = 1; + amrex::Print() << "Test failed. Max internal energy in cells too low: " << max_internal_energy << "\n"; + amrex::Print() << "Expected minimum max internal energy: " << expected_minimum_max_internal_energy << "\n"; + } + amrex::Print() << "---------------------------------------------------------" << "\n"; + + return status; +} diff --git a/src/problems/SN/test_SN.hpp b/src/problems/SN/test_SN.hpp new file mode 100644 index 0000000000..01baeb146f --- /dev/null +++ b/src/problems/SN/test_SN.hpp @@ -0,0 +1,10 @@ +#ifndef TEST_SN_HPP_ // NOLINT +#define TEST_SN_HPP_ +/// \file test_SN.hpp +/// \brief Defines a test problem for supernova feedback. +/// + +// function definitions +auto testproblem_SN() -> int; + +#endif // TEST_SN_HPP_ diff --git a/src/simulation.hpp b/src/simulation.hpp index 8845c3f55a..830fa72203 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -313,6 +313,8 @@ template class AMRSimulation : public amrex::AmrCore void incrementFluxRegisters(amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine, std::array &fluxArrays, int lev, amrex::Real dt_lev); + void particleMeshInteraction(amrex::Real time, amrex::Real dt); + // boundary condition AMREX_GPU_DEVICE static void setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &dest, int dcomp, int numcomp, amrex::GeometryData const &geom, amrex::Real time, const amrex::BCRec *bcr, int bcomp, @@ -1034,9 +1036,11 @@ template void AMRSimulation::evolve() // Only create particles at the finest level to avoid duplicate particle creation in regions where finer levels exist particleRegister_.createParticlesFromState(state_new_cc_[finest_level], finest_level, cur_time, dt_[0]); - // Stellar evolution and SN deposition - // TODO(cch): Need to take care of AMR subcycling - particleRegister_.depositSN(amrex::GetVecOfPtrs(state_new_cc_), 0, cur_time + dt_[0]); + // Stellar evolution and SN deposition; only apply to star particles + if (particleRegister_.HasStarParticles()) { + // TODO(cch): Need to take care of AMR subcycling + particleMeshInteraction(cur_time, dt_[0]); + } // Use the new type-aware particle destruction method // TODO(cch): Need to take care of AMR subcycling @@ -1382,6 +1386,23 @@ template void AMRSimulation::kickParticlesAllLev particleRegister_.kickParticlesAtLevel(lev, dt, accel_cc); } } + +template void AMRSimulation::particleMeshInteraction(amrex::Real time, amrex::Real dt) +{ + // Support up to 4 ghost cells for SN deposition. Default is 3, same as the TIGRESS model. + const int nghost = 3; + + // Assume all SN progenitors are at the finest level + const int lev = finest_level; + + // Create a MultiFab with 6 components (density, 3 x momentum, internal energy, energy) to hold the SN deposition + amrex::MultiFab state_buffer_at_level_cc(grids[lev], dmap[lev], Physics_NumVars::numHydroVars, nghost); + + // Deposit the SN particles into the MultiFab + particleRegister_.depositSN(state_new_cc_[lev], state_buffer_at_level_cc, lev, time, dt); + + // Sink accretion, to be implemented +} #endif // AMREX_SPACEDIM == 3 // N.B.: This function actually works for subcycled or not subcycled, as long as diff --git a/tests/SN.in b/tests/SN.in new file mode 100644 index 0000000000..fccc3a750a --- /dev/null +++ b/tests/SN.in @@ -0,0 +1,50 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = -3.949667304e+20 -3.949667304e+20 -3.949667304e+20 +geometry.prob_hi = 3.949667304e+20 3.949667304e+20 3.949667304e+20 +geometry.is_periodic = 1 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 1 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 32 32 32 +amr.max_level = 1 # number of levels = max_level + 1 +amr.blocking_factor = 32 # grid size must be divisible by this +amr.max_grid_size = 32 # at least 128 for GPUs +amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells +amr.grid_eff = 0.7 # default + +do_reflux = 1 +do_subcycle = 0 +do_tracers = 0 # turn on tracer particles + +hydro.artificial_viscosity_coefficient = 0.1 + +quokka.diagnostics = slice_z +quokka.slice_z.type = DiagFramePlane # Diagnostic type +quokka.slice_z.file = slicez_plt # Output file prefix (must end in "plt") +quokka.slice_z.normal = 2 # Plane normal (0 == x, 1 == y, 2 == z) +quokka.slice_z.center = 0.1 # Coordinate in the normal direction (cannot lie *exactly* on domain boundary) +quokka.slice_z.int = 1 # Output cadence (in number of coarse steps) +quokka.slice_z.interpolation = Linear # (Optional, default is Linear) Interpolation type: Linear or Quadratic +quokka.slice_z.field_names = gasDensity gasInternalEnergy x-GasMomentum y-GasMomentum z-GasMomentum # List of variables included in output + +cooling.enabled = 1 +cooling.cooling_table_type = "grackle" +cooling.hdf5_data_file = "../extern/grackle_data_files/input/CloudyData_UVB=HM2012.h5" +problem.SN_particles_file = "SN.txt" +problem.refine_half_domain = false + +problem.n_amb = 1.0 +problem.T_amb = 718.0 +particles.SN_scheme = 1 # 1: SN_thermal_or_thermal_momentum +particles.verbose = 1 + +max_timesteps = 10 +plotfile_interval = -1 diff --git a/tests/SN.txt b/tests/SN.txt new file mode 100644 index 0000000000..7e2b8652c9 --- /dev/null +++ b/tests/SN.txt @@ -0,0 +1,3 @@ +2 +0. 0. 0. 1.988409871e+34 0. 0. 0. 0. 0. 0. +3.949667304e+20 -3.949667304e+20 0. 1.988409871e+34 0. 0. 0. 0. 0. 0.