From 3b1b93ffe5dba7879bc60c781ac3a48b25a611b9 Mon Sep 17 00:00:00 2001 From: Piyush Sharda <34922596+psharda@users.noreply.github.com> Date: Fri, 19 Apr 2024 21:35:46 +0200 Subject: [PATCH 01/13] Print timestep info in scientific notation --- src/simulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation.hpp b/src/simulation.hpp index 200865f8c..a846d6ed3 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -1269,7 +1269,7 @@ template void AMRSimulation::timeStepWithSubcycl if (Verbose()) { amrex::Print() << "[Level " << lev << " step " << istep[lev] + 1 << "] "; - amrex::Print() << "ADVANCE with time = " << tNew_[lev] << " dt = " << dt_[lev] << '\n'; + amrex::Print() << "ADVANCE with time = " << std::scientific << tNew_[lev] << " dt = " << std::scientific << dt_[lev] << '\n'; } // Advance a single level for a single time step, and update flux registers From 09a713d1362aab21dab15dcab8e67c2bfe840fef Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 16:51:25 +0200 Subject: [PATCH 02/13] If burn fails, retry with reduced timestep --- src/Chemistry.hpp | 293 +++++++++++++++++++++++++--------------------- 1 file changed, 162 insertions(+), 131 deletions(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 47a396757..9360dc52e 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -26,146 +26,177 @@ namespace quokka::chemistry { +constexpr int max_retries = 5; // Maximum number of retries +constexpr Real dt_reduction_factor = 0.5; // Factor by which to reduce the timestep + AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt); -template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) +template void computeChemistry(amrex::MultiFab &mf, Real dt, const Real max_density_allowed, const Real min_density_allowed) { - - // Start off by assuming a successful burn. - int burn_success = 1; - - amrex::Gpu::Buffer d_num_failed({0}); - auto *p_num_failed = d_num_failed.data(); - - int num_failed = 0; - - const BL_PROFILE("Chemistry::computeChemistry()"); - for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { - const amrex::Box &indexRange = iter.validbox(); - auto const &state = mf.array(iter); - - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Real rho = state(i, j, k, HydroSystem::density_index); - const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); - const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); - const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); - const Real Ener = state(i, j, k, HydroSystem::energy_index); - const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); - - std::array chem = {-1.0}; - std::array inmfracs = {-1.0}; - Real insum = 0.0_rt; - - for (int nn = 0; nn < NumSpec; ++nn) { - chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / - rho; // state has partial densities, so divide by rho to get mass fractions - } - - // do chemistry using microphysics - - burn_t chemstate; - chemstate.success = true; - int burn_failed = 0; - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] = chem[nn] * rho / spmasses[nn]; - chemstate.xn[nn] = inmfracs[nn]; - } - - // dont do chemistry in cells with densities below the minimum density specified - if (rho < min_density_allowed) { - return; - } - - // stop the test if we have reached very high densities - if (rho > max_density_allowed) { - amrex::Abort("Density exceeded max_density_allowed!"); - } - - // input density and eint in burn state - // Microphysics needs specific eint - chemstate.rho = rho; - chemstate.e = Eint / rho; - - // call the EOS to set initial internal energy e - eos(eos_input_re, chemstate); - - // do the actual integration - // do it in .cpp so that it is not built at compile time for all tests - // which would otherwise slow down compilation due to the large RHS file - chemburner(chemstate, dt); - - if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { - amrex::Abort("Burner returned NAN"); - } - - if (!chemstate.success) { - burn_failed = 1; - } - - if (burn_failed) { - amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); - } - - // ensure positivity and normalize - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // update the number density of electrons due to charge conservation - // TODO(psharda): generalize this to other chem networks - chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + - chemstate.xn[9] + 2.0 * chemstate.xn[11]; - - // reconserve mass fractions post charge conservation - insum = 0; - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // get the updated specific eint - eos(eos_input_rt, chemstate); - - // get dEint - // Quokka uses rho*eint - const Real dEint = (chemstate.e * chemstate.rho) - Eint; - state(i, j, k, HydroSystem::internalEnergy_index) += dEint; - state(i, j, k, HydroSystem::energy_index) += dEint; - - for (int nn = 0; nn < NumSpec; ++nn) { - state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities - } - }); + // Start off by assuming a successful burn. + int burn_success = 1; + + amrex::Gpu::Buffer d_num_failed({0}); + auto *p_num_failed = d_num_failed.data(); + + int num_failed = 0; + + const BL_PROFILE("Chemistry::computeChemistry()"); + + // Loop for retries + for (int retry = 0; retry < max_retries; ++retry) + { + burn_success = 1; // Reset burn_success flag + + // Loop over MultiFab + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) + { + const amrex::Box &indexRange = iter.validbox(); + auto const &state = mf.array(iter); + + // Perform burn computation + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const Real rho = state(i, j, k, HydroSystem::density_index); + const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); + const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); + const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); + const Real Ener = state(i, j, k, HydroSystem::energy_index); + const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); + + std::array chem = {-1.0}; + std::array inmfracs = {-1.0}; + Real insum = 0.0_rt; + + for (int nn = 0; nn < NumSpec; ++nn) { + chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / + rho; // state has partial densities, so divide by rho to get mass fractions + } + + // do chemistry using microphysics + + burn_t chemstate; + chemstate.success = true; + int burn_failed = 0; + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] = chem[nn] * rho / spmasses[nn]; + chemstate.xn[nn] = inmfracs[nn]; + } + + // dont do chemistry in cells with densities below the minimum density specified + if (rho < min_density_allowed) { + return; + } + + // stop the test if we have reached very high densities + if (rho > max_density_allowed) { + amrex::Abort("Density exceeded max_density_allowed!"); + } + + // input density and eint in burn state + // Microphysics needs specific eint + chemstate.rho = rho; + chemstate.e = Eint / rho; + + // call the EOS to set initial internal energy e + eos(eos_input_re, chemstate); + + // do the actual integration + // do it in .cpp so that it is not built at compile time for all tests + // which would otherwise slow down compilation due to the large RHS file + chemburner(chemstate, dt); + + if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { + amrex::Abort("Burner returned NAN"); + } + + if (!chemstate.success) { + burn_failed = 1; + } + + if (burn_failed) { + amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); + } + + // ensure positivity and normalize + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // update the number density of electrons due to charge conservation + // TODO(psharda): generalize this to other chem networks + chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + + chemstate.xn[9] + 2.0 * chemstate.xn[11]; + + // reconserve mass fractions post charge conservation + insum = 0; + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // get the updated specific eint + eos(eos_input_rt, chemstate); + + // get dEint + // Quokka uses rho*eint + const Real dEint = (chemstate.e * chemstate.rho) - Eint; + state(i, j, k, HydroSystem::internalEnergy_index) += dEint; + state(i, j, k, HydroSystem::energy_index) += dEint; + + for (int nn = 0; nn < NumSpec; ++nn) { + state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities + } + }); #if defined(AMREX_USE_HIP) - amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. + amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. #endif - } + } + + num_failed = *(d_num_failed.copyToHost()); + + burn_success = !num_failed; + amrex::ParallelDescriptor::ReduceIntMin(burn_success); + + // If burn was successful, break out of the retry loop + if (burn_success) + break; + + // If burn was not successful and we haven't reached max retries, reduce the timestep + if (retry < max_retries - 1) + { + // Reduce timestep + Real new_dt = dt * dt_reduction_factor; - num_failed = *(d_num_failed.copyToHost()); + // Output information about reduced timestep + amrex::Print() << "Burn failed. Reducing timestep to " << new_dt << std::endl; - burn_success = !num_failed; - amrex::ParallelDescriptor::ReduceIntMin(burn_success); + // Update dt for the next iteration + dt = new_dt; + } + } - if (!burn_success) { - amrex::Abort("Burn failed in VODE. Aborting."); - } + // If burn was not successful after max retries, abort + if (!burn_success) + { + amrex::Abort("Burn failed after max retries."); + } } } // namespace quokka::chemistry From 92bc05329cd27bdd14299c4dc15b693735a33617 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 20 Apr 2024 14:55:31 +0000 Subject: [PATCH 03/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/Chemistry.hpp | 307 +++++++++++++++++++++++----------------------- 1 file changed, 152 insertions(+), 155 deletions(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 9360dc52e..251261c9c 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -26,177 +26,174 @@ namespace quokka::chemistry { -constexpr int max_retries = 5; // Maximum number of retries +constexpr int max_retries = 5; // Maximum number of retries constexpr Real dt_reduction_factor = 0.5; // Factor by which to reduce the timestep AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt); template void computeChemistry(amrex::MultiFab &mf, Real dt, const Real max_density_allowed, const Real min_density_allowed) { - // Start off by assuming a successful burn. - int burn_success = 1; - - amrex::Gpu::Buffer d_num_failed({0}); - auto *p_num_failed = d_num_failed.data(); - - int num_failed = 0; - - const BL_PROFILE("Chemistry::computeChemistry()"); - - // Loop for retries - for (int retry = 0; retry < max_retries; ++retry) - { - burn_success = 1; // Reset burn_success flag - - // Loop over MultiFab - for (amrex::MFIter iter(mf); iter.isValid(); ++iter) - { - const amrex::Box &indexRange = iter.validbox(); - auto const &state = mf.array(iter); - - // Perform burn computation - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Real rho = state(i, j, k, HydroSystem::density_index); - const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); - const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); - const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); - const Real Ener = state(i, j, k, HydroSystem::energy_index); - const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); - - std::array chem = {-1.0}; - std::array inmfracs = {-1.0}; - Real insum = 0.0_rt; - - for (int nn = 0; nn < NumSpec; ++nn) { - chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / - rho; // state has partial densities, so divide by rho to get mass fractions - } - - // do chemistry using microphysics - - burn_t chemstate; - chemstate.success = true; - int burn_failed = 0; - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] = chem[nn] * rho / spmasses[nn]; - chemstate.xn[nn] = inmfracs[nn]; - } - - // dont do chemistry in cells with densities below the minimum density specified - if (rho < min_density_allowed) { - return; - } - - // stop the test if we have reached very high densities - if (rho > max_density_allowed) { - amrex::Abort("Density exceeded max_density_allowed!"); - } - - // input density and eint in burn state - // Microphysics needs specific eint - chemstate.rho = rho; - chemstate.e = Eint / rho; - - // call the EOS to set initial internal energy e - eos(eos_input_re, chemstate); - - // do the actual integration - // do it in .cpp so that it is not built at compile time for all tests - // which would otherwise slow down compilation due to the large RHS file - chemburner(chemstate, dt); - - if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { - amrex::Abort("Burner returned NAN"); - } - - if (!chemstate.success) { - burn_failed = 1; - } - - if (burn_failed) { - amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); - } - - // ensure positivity and normalize - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // update the number density of electrons due to charge conservation - // TODO(psharda): generalize this to other chem networks - chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + - chemstate.xn[9] + 2.0 * chemstate.xn[11]; - - // reconserve mass fractions post charge conservation - insum = 0; - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // get the updated specific eint - eos(eos_input_rt, chemstate); - - // get dEint - // Quokka uses rho*eint - const Real dEint = (chemstate.e * chemstate.rho) - Eint; - state(i, j, k, HydroSystem::internalEnergy_index) += dEint; - state(i, j, k, HydroSystem::energy_index) += dEint; - - for (int nn = 0; nn < NumSpec; ++nn) { - state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities - } - }); + // Start off by assuming a successful burn. + int burn_success = 1; + + amrex::Gpu::Buffer d_num_failed({0}); + auto *p_num_failed = d_num_failed.data(); + + int num_failed = 0; + + const BL_PROFILE("Chemistry::computeChemistry()"); + + // Loop for retries + for (int retry = 0; retry < max_retries; ++retry) { + burn_success = 1; // Reset burn_success flag + + // Loop over MultiFab + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); + auto const &state = mf.array(iter); + + // Perform burn computation + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const Real rho = state(i, j, k, HydroSystem::density_index); + const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); + const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); + const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); + const Real Ener = state(i, j, k, HydroSystem::energy_index); + const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); + + std::array chem = {-1.0}; + std::array inmfracs = {-1.0}; + Real insum = 0.0_rt; + + for (int nn = 0; nn < NumSpec; ++nn) { + chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / + rho; // state has partial densities, so divide by rho to get mass fractions + } + + // do chemistry using microphysics + + burn_t chemstate; + chemstate.success = true; + int burn_failed = 0; + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] = chem[nn] * rho / spmasses[nn]; + chemstate.xn[nn] = inmfracs[nn]; + } + + // dont do chemistry in cells with densities below the minimum density specified + if (rho < min_density_allowed) { + return; + } + + // stop the test if we have reached very high densities + if (rho > max_density_allowed) { + amrex::Abort("Density exceeded max_density_allowed!"); + } + + // input density and eint in burn state + // Microphysics needs specific eint + chemstate.rho = rho; + chemstate.e = Eint / rho; + + // call the EOS to set initial internal energy e + eos(eos_input_re, chemstate); + + // do the actual integration + // do it in .cpp so that it is not built at compile time for all tests + // which would otherwise slow down compilation due to the large RHS file + chemburner(chemstate, dt); + + if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { + amrex::Abort("Burner returned NAN"); + } + + if (!chemstate.success) { + burn_failed = 1; + } + + if (burn_failed) { + amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); + } + + // ensure positivity and normalize + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // update the number density of electrons due to charge conservation + // TODO(psharda): generalize this to other chem networks + chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + + chemstate.xn[9] + 2.0 * chemstate.xn[11]; + + // reconserve mass fractions post charge conservation + insum = 0; + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // get the updated specific eint + eos(eos_input_rt, chemstate); + + // get dEint + // Quokka uses rho*eint + const Real dEint = (chemstate.e * chemstate.rho) - Eint; + state(i, j, k, HydroSystem::internalEnergy_index) += dEint; + state(i, j, k, HydroSystem::energy_index) += dEint; + + for (int nn = 0; nn < NumSpec; ++nn) { + state(i, j, k, HydroSystem::scalar0_index + nn) = + inmfracs[nn] * rho; // scale by rho to return partial densities + } + }); #if defined(AMREX_USE_HIP) - amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. + amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. #endif - } + } - num_failed = *(d_num_failed.copyToHost()); + num_failed = *(d_num_failed.copyToHost()); - burn_success = !num_failed; - amrex::ParallelDescriptor::ReduceIntMin(burn_success); + burn_success = !num_failed; + amrex::ParallelDescriptor::ReduceIntMin(burn_success); - // If burn was successful, break out of the retry loop - if (burn_success) - break; + // If burn was successful, break out of the retry loop + if (burn_success) + break; - // If burn was not successful and we haven't reached max retries, reduce the timestep - if (retry < max_retries - 1) - { - // Reduce timestep - Real new_dt = dt * dt_reduction_factor; + // If burn was not successful and we haven't reached max retries, reduce the timestep + if (retry < max_retries - 1) { + // Reduce timestep + Real new_dt = dt * dt_reduction_factor; - // Output information about reduced timestep - amrex::Print() << "Burn failed. Reducing timestep to " << new_dt << std::endl; + // Output information about reduced timestep + amrex::Print() << "Burn failed. Reducing timestep to " << new_dt << std::endl; - // Update dt for the next iteration - dt = new_dt; - } - } + // Update dt for the next iteration + dt = new_dt; + } + } - // If burn was not successful after max retries, abort - if (!burn_success) - { - amrex::Abort("Burn failed after max retries."); - } + // If burn was not successful after max retries, abort + if (!burn_success) { + amrex::Abort("Burn failed after max retries."); + } } } // namespace quokka::chemistry From 0a78ca9ab66e4944c5af4747f32b22f059827aa4 Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 17:00:12 +0200 Subject: [PATCH 04/13] add braces --- src/Chemistry.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 251261c9c..30383b27d 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -174,8 +174,9 @@ template void computeChemistry(amrex::MultiFab &mf, Real dt amrex::ParallelDescriptor::ReduceIntMin(burn_success); // If burn was successful, break out of the retry loop - if (burn_success) + if (burn_success) { break; + } // If burn was not successful and we haven't reached max retries, reduce the timestep if (retry < max_retries - 1) { From 11e906bfabe0090d191de45ac98743fe98977e1d Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 21:24:47 +0200 Subject: [PATCH 05/13] reset file --- src/Chemistry.hpp | 291 +++++++++++++++++++++------------------------- 1 file changed, 131 insertions(+), 160 deletions(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 30383b27d..2f6544867 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -26,175 +26,146 @@ namespace quokka::chemistry { -constexpr int max_retries = 5; // Maximum number of retries -constexpr Real dt_reduction_factor = 0.5; // Factor by which to reduce the timestep - AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt); -template void computeChemistry(amrex::MultiFab &mf, Real dt, const Real max_density_allowed, const Real min_density_allowed) +template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) { - // Start off by assuming a successful burn. - int burn_success = 1; - - amrex::Gpu::Buffer d_num_failed({0}); - auto *p_num_failed = d_num_failed.data(); - - int num_failed = 0; - - const BL_PROFILE("Chemistry::computeChemistry()"); - - // Loop for retries - for (int retry = 0; retry < max_retries; ++retry) { - burn_success = 1; // Reset burn_success flag - - // Loop over MultiFab - for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { - const amrex::Box &indexRange = iter.validbox(); - auto const &state = mf.array(iter); - - // Perform burn computation - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Real rho = state(i, j, k, HydroSystem::density_index); - const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); - const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); - const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); - const Real Ener = state(i, j, k, HydroSystem::energy_index); - const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); - - std::array chem = {-1.0}; - std::array inmfracs = {-1.0}; - Real insum = 0.0_rt; - - for (int nn = 0; nn < NumSpec; ++nn) { - chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / - rho; // state has partial densities, so divide by rho to get mass fractions - } - - // do chemistry using microphysics - - burn_t chemstate; - chemstate.success = true; - int burn_failed = 0; - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] = chem[nn] * rho / spmasses[nn]; - chemstate.xn[nn] = inmfracs[nn]; - } - - // dont do chemistry in cells with densities below the minimum density specified - if (rho < min_density_allowed) { - return; - } - - // stop the test if we have reached very high densities - if (rho > max_density_allowed) { - amrex::Abort("Density exceeded max_density_allowed!"); - } - - // input density and eint in burn state - // Microphysics needs specific eint - chemstate.rho = rho; - chemstate.e = Eint / rho; - - // call the EOS to set initial internal energy e - eos(eos_input_re, chemstate); - - // do the actual integration - // do it in .cpp so that it is not built at compile time for all tests - // which would otherwise slow down compilation due to the large RHS file - chemburner(chemstate, dt); - - if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { - amrex::Abort("Burner returned NAN"); - } - - if (!chemstate.success) { - burn_failed = 1; - } - - if (burn_failed) { - amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); - } - - // ensure positivity and normalize - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // update the number density of electrons due to charge conservation - // TODO(psharda): generalize this to other chem networks - chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + - chemstate.xn[9] + 2.0 * chemstate.xn[11]; - - // reconserve mass fractions post charge conservation - insum = 0; - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // get the updated specific eint - eos(eos_input_rt, chemstate); - - // get dEint - // Quokka uses rho*eint - const Real dEint = (chemstate.e * chemstate.rho) - Eint; - state(i, j, k, HydroSystem::internalEnergy_index) += dEint; - state(i, j, k, HydroSystem::energy_index) += dEint; - - for (int nn = 0; nn < NumSpec; ++nn) { - state(i, j, k, HydroSystem::scalar0_index + nn) = - inmfracs[nn] * rho; // scale by rho to return partial densities - } - }); + + // Start off by assuming a successful burn. + int burn_success = 1; + + amrex::Gpu::Buffer d_num_failed({0}); + auto *p_num_failed = d_num_failed.data(); + + int num_failed = 0; + + const BL_PROFILE("Chemistry::computeChemistry()"); + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); + auto const &state = mf.array(iter); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const Real rho = state(i, j, k, HydroSystem::density_index); + const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); + const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); + const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); + const Real Ener = state(i, j, k, HydroSystem::energy_index); + const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); + + std::array chem = {-1.0}; + std::array inmfracs = {-1.0}; + Real insum = 0.0_rt; + + for (int nn = 0; nn < NumSpec; ++nn) { + chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / + rho; // state has partial densities, so divide by rho to get mass fractions + } + + // do chemistry using microphysics + + burn_t chemstate; + chemstate.success = true; + int burn_failed = 0; + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] = chem[nn] * rho / spmasses[nn]; + chemstate.xn[nn] = inmfracs[nn]; + } + + // dont do chemistry in cells with densities below the minimum density specified + if (rho < min_density_allowed) { + return; + } + + // stop the test if we have reached very high densities + if (rho > max_density_allowed) { + amrex::Abort("Density exceeded max_density_allowed!"); + } + + // input density and eint in burn state + // Microphysics needs specific eint + chemstate.rho = rho; + chemstate.e = Eint / rho; + + // call the EOS to set initial internal energy e + eos(eos_input_re, chemstate); + + // do the actual integration + // do it in .cpp so that it is not built at compile time for all tests + // which would otherwise slow down compilation due to the large RHS file + chemburner(chemstate, dt); + + if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { + amrex::Abort("Burner returned NAN"); + } + + if (!chemstate.success) { + burn_failed = 1; + } + + if (burn_failed) { + amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); + } + + // ensure positivity and normalize + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // update the number density of electrons due to charge conservation + // TODO(psharda): generalize this to other chem networks + chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + + chemstate.xn[9] + 2.0 * chemstate.xn[11]; + + // reconserve mass fractions post charge conservation + insum = 0; + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // get the updated specific eint + eos(eos_input_rt, chemstate); + + // get dEint + // Quokka uses rho*eint + const Real dEint = (chemstate.e * chemstate.rho) - Eint; + state(i, j, k, HydroSystem::internalEnergy_index) += dEint; + state(i, j, k, HydroSystem::energy_index) += dEint; + + for (int nn = 0; nn < NumSpec; ++nn) { + state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities + } + }); #if defined(AMREX_USE_HIP) - amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. + amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. #endif - } - - num_failed = *(d_num_failed.copyToHost()); - - burn_success = !num_failed; - amrex::ParallelDescriptor::ReduceIntMin(burn_success); - - // If burn was successful, break out of the retry loop - if (burn_success) { - break; - } - - // If burn was not successful and we haven't reached max retries, reduce the timestep - if (retry < max_retries - 1) { - // Reduce timestep - Real new_dt = dt * dt_reduction_factor; + } - // Output information about reduced timestep - amrex::Print() << "Burn failed. Reducing timestep to " << new_dt << std::endl; + num_failed = *(d_num_failed.copyToHost()); - // Update dt for the next iteration - dt = new_dt; - } - } + burn_success = !num_failed; + amrex::ParallelDescriptor::ReduceIntMin(burn_success); - // If burn was not successful after max retries, abort - if (!burn_success) { - amrex::Abort("Burn failed after max retries."); - } + if (!burn_success) { + amrex::Abort("Burn failed in VODE. Aborting."); + } } } // namespace quokka::chemistry From 527f15e62448d53a13e4dd89c42db108a1a1a8f3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 20 Apr 2024 19:25:05 +0000 Subject: [PATCH 06/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/Chemistry.hpp | 258 +++++++++++++++++++++++----------------------- 1 file changed, 129 insertions(+), 129 deletions(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 2f6544867..47a396757 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -31,141 +31,141 @@ AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt); template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) { - // Start off by assuming a successful burn. - int burn_success = 1; - - amrex::Gpu::Buffer d_num_failed({0}); - auto *p_num_failed = d_num_failed.data(); - - int num_failed = 0; - - const BL_PROFILE("Chemistry::computeChemistry()"); - for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { - const amrex::Box &indexRange = iter.validbox(); - auto const &state = mf.array(iter); - - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Real rho = state(i, j, k, HydroSystem::density_index); - const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); - const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); - const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); - const Real Ener = state(i, j, k, HydroSystem::energy_index); - const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); - - std::array chem = {-1.0}; - std::array inmfracs = {-1.0}; - Real insum = 0.0_rt; - - for (int nn = 0; nn < NumSpec; ++nn) { - chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / - rho; // state has partial densities, so divide by rho to get mass fractions - } - - // do chemistry using microphysics - - burn_t chemstate; - chemstate.success = true; - int burn_failed = 0; - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] = chem[nn] * rho / spmasses[nn]; - chemstate.xn[nn] = inmfracs[nn]; - } - - // dont do chemistry in cells with densities below the minimum density specified - if (rho < min_density_allowed) { - return; - } - - // stop the test if we have reached very high densities - if (rho > max_density_allowed) { - amrex::Abort("Density exceeded max_density_allowed!"); - } - - // input density and eint in burn state - // Microphysics needs specific eint - chemstate.rho = rho; - chemstate.e = Eint / rho; - - // call the EOS to set initial internal energy e - eos(eos_input_re, chemstate); - - // do the actual integration - // do it in .cpp so that it is not built at compile time for all tests - // which would otherwise slow down compilation due to the large RHS file - chemburner(chemstate, dt); - - if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { - amrex::Abort("Burner returned NAN"); - } - - if (!chemstate.success) { - burn_failed = 1; - } - - if (burn_failed) { - amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); - } - - // ensure positivity and normalize - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // update the number density of electrons due to charge conservation - // TODO(psharda): generalize this to other chem networks - chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + - chemstate.xn[9] + 2.0 * chemstate.xn[11]; - - // reconserve mass fractions post charge conservation - insum = 0; - for (int nn = 0; nn < NumSpec; ++nn) { - chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); - inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; - insum += inmfracs[nn]; - } - - for (int nn = 0; nn < NumSpec; ++nn) { - inmfracs[nn] /= insum; - // update the number densities with conserved mass fractions - chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; - } - - // get the updated specific eint - eos(eos_input_rt, chemstate); - - // get dEint - // Quokka uses rho*eint - const Real dEint = (chemstate.e * chemstate.rho) - Eint; - state(i, j, k, HydroSystem::internalEnergy_index) += dEint; - state(i, j, k, HydroSystem::energy_index) += dEint; - - for (int nn = 0; nn < NumSpec; ++nn) { - state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities - } - }); + // Start off by assuming a successful burn. + int burn_success = 1; + + amrex::Gpu::Buffer d_num_failed({0}); + auto *p_num_failed = d_num_failed.data(); + + int num_failed = 0; + + const BL_PROFILE("Chemistry::computeChemistry()"); + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); + auto const &state = mf.array(iter); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const Real rho = state(i, j, k, HydroSystem::density_index); + const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); + const Real ymom = state(i, j, k, HydroSystem::x2Momentum_index); + const Real zmom = state(i, j, k, HydroSystem::x3Momentum_index); + const Real Ener = state(i, j, k, HydroSystem::energy_index); + const Real Eint = RadSystem::ComputeEintFromEgas(rho, xmom, ymom, zmom, Ener); + + std::array chem = {-1.0}; + std::array inmfracs = {-1.0}; + Real insum = 0.0_rt; + + for (int nn = 0; nn < NumSpec; ++nn) { + chem[nn] = state(i, j, k, HydroSystem::scalar0_index + nn) / + rho; // state has partial densities, so divide by rho to get mass fractions + } + + // do chemistry using microphysics + + burn_t chemstate; + chemstate.success = true; + int burn_failed = 0; + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] = chem[nn] * rho / spmasses[nn]; + chemstate.xn[nn] = inmfracs[nn]; + } + + // dont do chemistry in cells with densities below the minimum density specified + if (rho < min_density_allowed) { + return; + } + + // stop the test if we have reached very high densities + if (rho > max_density_allowed) { + amrex::Abort("Density exceeded max_density_allowed!"); + } + + // input density and eint in burn state + // Microphysics needs specific eint + chemstate.rho = rho; + chemstate.e = Eint / rho; + + // call the EOS to set initial internal energy e + eos(eos_input_re, chemstate); + + // do the actual integration + // do it in .cpp so that it is not built at compile time for all tests + // which would otherwise slow down compilation due to the large RHS file + chemburner(chemstate, dt); + + if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) { + amrex::Abort("Burner returned NAN"); + } + + if (!chemstate.success) { + burn_failed = 1; + } + + if (burn_failed) { + amrex::Gpu::Atomic::Add(p_num_failed, burn_failed); + } + + // ensure positivity and normalize + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // update the number density of electrons due to charge conservation + // TODO(psharda): generalize this to other chem networks + chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] + + chemstate.xn[9] + 2.0 * chemstate.xn[11]; + + // reconserve mass fractions post charge conservation + insum = 0; + for (int nn = 0; nn < NumSpec; ++nn) { + chemstate.xn[nn] = amrex::max(chemstate.xn[nn], small_x); + inmfracs[nn] = spmasses[nn] * chemstate.xn[nn] / chemstate.rho; + insum += inmfracs[nn]; + } + + for (int nn = 0; nn < NumSpec; ++nn) { + inmfracs[nn] /= insum; + // update the number densities with conserved mass fractions + chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn]; + } + + // get the updated specific eint + eos(eos_input_rt, chemstate); + + // get dEint + // Quokka uses rho*eint + const Real dEint = (chemstate.e * chemstate.rho) - Eint; + state(i, j, k, HydroSystem::internalEnergy_index) += dEint; + state(i, j, k, HydroSystem::energy_index) += dEint; + + for (int nn = 0; nn < NumSpec; ++nn) { + state(i, j, k, HydroSystem::scalar0_index + nn) = inmfracs[nn] * rho; // scale by rho to return partial densities + } + }); #if defined(AMREX_USE_HIP) - amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. + amrex::Gpu::streamSynchronize(); // otherwise HIP may fail to allocate the necessary resources. #endif - } + } - num_failed = *(d_num_failed.copyToHost()); + num_failed = *(d_num_failed.copyToHost()); - burn_success = !num_failed; - amrex::ParallelDescriptor::ReduceIntMin(burn_success); + burn_success = !num_failed; + amrex::ParallelDescriptor::ReduceIntMin(burn_success); - if (!burn_success) { - amrex::Abort("Burn failed in VODE. Aborting."); - } + if (!burn_success) { + amrex::Abort("Burn failed in VODE. Aborting."); + } } } // namespace quokka::chemistry From 7888ac4f1531e1960fa2c487138497ed785cd7a1 Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 22:08:02 +0200 Subject: [PATCH 07/13] return burn_success as a boolean --- src/Chemistry.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 47a396757..50392553f 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -28,7 +28,7 @@ namespace quokka::chemistry AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt); -template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) +template auto computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) -> bool { // Start off by assuming a successful burn. @@ -164,8 +164,11 @@ template void computeChemistry(amrex::MultiFab &mf, const R amrex::ParallelDescriptor::ReduceIntMin(burn_success); if (!burn_success) { - amrex::Abort("Burn failed in VODE. Aborting."); + // amrex::Abort("Burn failed in VODE. Aborting."); + amrex::Print() << "WARNNING: Unsuccessful burn. Retrying hydro step." << "\n"; } + + return burn_success; } } // namespace quokka::chemistry From 72aaa8fe714f5bbd6d5287fe0ab633d73234b83c Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 22:08:24 +0200 Subject: [PATCH 08/13] retry hydro if chemistry burn failed --- src/RadhydroSimulation.hpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index c1b5654ed..8bbd38eb2 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -230,7 +230,7 @@ template class RadhydroSimulation : public AMRSimulation bool; void addStrangSplitSources(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt_lev); - void addStrangSplitSourcesWithBuiltin(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt_lev); + auto addStrangSplitSourcesWithBuiltin(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt_lev) -> bool; auto isCflViolated(int lev, amrex::Real time, amrex::Real dt_actual) -> bool; @@ -507,8 +507,12 @@ template void RadhydroSimulation::addStrangSplit } template -void RadhydroSimulation::addStrangSplitSourcesWithBuiltin(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt) +auto RadhydroSimulation::addStrangSplitSourcesWithBuiltin(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt) -> bool { + + // start by assuming chemistry burn is successful. + bool burn_success = true; + if (enableCooling_ == 1) { // compute cooling if (coolingTableType_ == "grackle") { @@ -523,12 +527,14 @@ void RadhydroSimulation::addStrangSplitSourcesWithBuiltin(amrex::Mult #ifdef PRIMORDIAL_CHEM if (enableChemistry_ == 1) { // compute chemistry - quokka::chemistry::computeChemistry(state, dt, max_density_allowed, min_density_allowed); + burn_success = quokka::chemistry::computeChemistry(state, dt, max_density_allowed, min_density_allowed); } #endif // compute user-specified sources addStrangSplitSources(state, lev, time, dt); + + return burn_success; } template @@ -1028,7 +1034,7 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o auto dx = geom[lev].CellSizeArray(); // do Strang split source terms (first half-step) - addStrangSplitSourcesWithBuiltin(state_old_cc_tmp, lev, time, 0.5 * dt_lev); + auto burn_success_first = addStrangSplitSourcesWithBuiltin(state_old_cc_tmp, lev, time, 0.5 * dt_lev); // create temporary multifab for intermediate state amrex::MultiFab state_inter_cc_(grids[lev], dmap[lev], Physics_Indices::nvarTotal_cc, nghost_cc_); @@ -1290,10 +1296,10 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o #endif // do Strang split source terms (second half-step) - addStrangSplitSourcesWithBuiltin(state_new_cc_[lev], lev, time + dt_lev, 0.5 * dt_lev); + auto burn_success_second = addStrangSplitSourcesWithBuiltin(state_new_cc_[lev], lev, time + dt_lev, 0.5 * dt_lev); - // check if we have violated the CFL timestep - return !isCflViolated(lev, time, dt_lev); + // check if we have violated the CFL timestep or burn failed in chemistry + return (!isCflViolated(lev, time, dt_lev) && burn_success_first && burn_success_second); } template From 41bbe803402502299479d34a19927465fe4efc86 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 20 Apr 2024 20:08:39 +0000 Subject: [PATCH 09/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/Chemistry.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 50392553f..b387d2a9f 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -165,7 +165,8 @@ template auto computeChemistry(amrex::MultiFab &mf, const R if (!burn_success) { // amrex::Abort("Burn failed in VODE. Aborting."); - amrex::Print() << "WARNNING: Unsuccessful burn. Retrying hydro step." << "\n"; + amrex::Print() << "WARNNING: Unsuccessful burn. Retrying hydro step." + << "\n"; } return burn_success; From 7d0c3f52e56e3b60558235d22ff776cb6accc127 Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 22:38:47 +0200 Subject: [PATCH 10/13] silence clang-tidy --- src/RadhydroSimulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 8bbd38eb2..7d13f73dc 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -511,7 +511,7 @@ auto RadhydroSimulation::addStrangSplitSourcesWithBuiltin(amrex::Mult { // start by assuming chemistry burn is successful. - bool burn_success = true; + bool burn_success = true; // NOLINT if (enableCooling_ == 1) { // compute cooling From 0eae594dc459d80fd25c361a91e81bf307c6a8ac Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sat, 20 Apr 2024 23:07:00 +0200 Subject: [PATCH 11/13] correct spelling --- src/Chemistry.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index b387d2a9f..0e289148c 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -165,7 +165,7 @@ template auto computeChemistry(amrex::MultiFab &mf, const R if (!burn_success) { // amrex::Abort("Burn failed in VODE. Aborting."); - amrex::Print() << "WARNNING: Unsuccessful burn. Retrying hydro step." + amrex::Print() << "\t>> WARNING: Unsuccessful burn. Retrying hydro step." << "\n"; } From 45e091886c52c4b9bcbc13e7b300f663f9fb21bd Mon Sep 17 00:00:00 2001 From: Piyush Sharda Date: Sun, 21 Apr 2024 17:03:12 +0200 Subject: [PATCH 12/13] exit already if first burn failed --- src/RadhydroSimulation.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 7d13f73dc..1186f99f3 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -1036,6 +1036,12 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o // do Strang split source terms (first half-step) auto burn_success_first = addStrangSplitSourcesWithBuiltin(state_old_cc_tmp, lev, time, 0.5 * dt_lev); + // check if burn failed in chemistry. If it did, return + if (!burn_success_first) { + return burn_success_first; + } + + // create temporary multifab for intermediate state amrex::MultiFab state_inter_cc_(grids[lev], dmap[lev], Physics_Indices::nvarTotal_cc, nghost_cc_); state_inter_cc_.setVal(0); // prevent assert in fillBoundaryConditions when radiation is enabled @@ -1299,7 +1305,7 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o auto burn_success_second = addStrangSplitSourcesWithBuiltin(state_new_cc_[lev], lev, time + dt_lev, 0.5 * dt_lev); // check if we have violated the CFL timestep or burn failed in chemistry - return (!isCflViolated(lev, time, dt_lev) && burn_success_first && burn_success_second); + return (!isCflViolated(lev, time, dt_lev) && burn_success_second); } template From 8f8380446989abfadc0ab1031b767aa78cc057cf Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 21 Apr 2024 15:03:36 +0000 Subject: [PATCH 13/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/RadhydroSimulation.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 1186f99f3..92345d849 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -1041,7 +1041,6 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o return burn_success_first; } - // create temporary multifab for intermediate state amrex::MultiFab state_inter_cc_(grids[lev], dmap[lev], Physics_Indices::nvarTotal_cc, nghost_cc_); state_inter_cc_.setVal(0); // prevent assert in fillBoundaryConditions when radiation is enabled