From 86686408835cf731c5792151bb890f6c6f78b2fc Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Tue, 20 Aug 2024 02:04:57 +1000 Subject: [PATCH] Make assertion work in the exchange step and dust temperature solution step (#716) ### Description A fix to #543 : make the Newton-Raphson iteration return a flag and check it on host. More changes: - Change `maxIter` from 400 to 50. 50 iterations should be more than enough. If it doesn't converge in less than 50 steps, it probably never will. ### Related issues Fixes #543 ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- src/QuokkaSimulation.hpp | 38 ++++++++++++++++++++++++++---- src/radiation/radiation_system.hpp | 36 +++++++++++++++++++++------- 2 files changed, 60 insertions(+), 14 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index a1278dc44..12218ae1d 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -251,7 +251,8 @@ template class QuokkaSimulation : public AMRSimulation const &stateNew, const amrex::Box &indexRange, amrex::Real time, double dt, int stage, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, - amrex::GpuArray const &prob_hi); + amrex::GpuArray const &prob_hi, int *num_failed_coupling, int *num_failed_dust, + int *p_num_failed_outer_ite); auto computeRadiationFluxes(amrex::Array4 const &consVar, const amrex::Box &indexRange, int nvars, amrex::GpuArray dx) @@ -1615,8 +1616,16 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real // new radiation state is stored in state_new_cc_ // new hydro state is stored in state_new_cc_ (always the case during radiation update) + amrex::Gpu::Buffer num_failed_coupling({0}); + amrex::Gpu::Buffer num_failed_dust({0}); + amrex::Gpu::Buffer num_failed_outer({0}); + int *p_num_failed_coupling = num_failed_coupling.data(); + int *p_num_failed_dust = num_failed_dust.data(); + int *p_num_failed_outer = num_failed_outer.data(); + if constexpr (IMEX_a22 > 0.0) { // matter-radiation exchange source terms of stage 1 + for (amrex::MFIter iter(state_new_cc_[lev]); iter.isValid(); ++iter) { const amrex::Box &indexRange = iter.validbox(); auto const &stateNew = state_new_cc_[lev].array(iter); @@ -1625,7 +1634,8 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real // update state_new_cc_[lev] in place (updates both radiation and hydro vars) // Note that only a fraction (IMEX_a32) of the matter-radiation exchange source terms are added to hydro. This ensures that the // hydro properties get to t + IMEX_a32 dt in terms of matter-radiation exchange. - operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi); + operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi, p_num_failed_coupling, + p_num_failed_dust, p_num_failed_outer); } } @@ -1642,7 +1652,23 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real auto const &prob_lo = geom[lev].ProbLoArray(); auto const &prob_hi = geom[lev].ProbHiArray(); // update state_new_cc_[lev] in place (updates both radiation and hydro vars) - operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi); + operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi, p_num_failed_coupling, + p_num_failed_dust, p_num_failed_outer); + } + + const int nf_coupling = *(num_failed_coupling.copyToHost()); + const int nf_dust = *(num_failed_dust.copyToHost()); + const int nf_outer = *(num_failed_outer.copyToHost()); + // Note that the nf_dust has to abort BEFORE nf_coupling, because the dust temperature is used in the matter-radiation coupling and if dust + // temperature is negative, the matter-radiation coupling will fail to converge. + if (nf_dust > 0) { + amrex::Abort("Newton-Raphson iteration for dust temperature failed to converge or dust temperature is negative!"); + } + if (nf_coupling > 0) { + amrex::Abort("Newton-Raphson iteration for matter-radiation coupling failed to converge!"); + } + if (nf_outer > 0) { + amrex::Abort("Outer iteration for matter-radiation coupling failed to converge!"); } // new hydro+radiation state is stored in state_new_cc_ @@ -1794,7 +1820,8 @@ template void QuokkaSimulation::operatorSplitSourceTerms(amrex::Array4 const &stateNew, const amrex::Box &indexRange, const amrex::Real time, const double dt, const int stage, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, - amrex::GpuArray const &prob_hi) + amrex::GpuArray const &prob_hi, int *p_num_failed_coupling, + int *p_num_failed_dust, int *p_num_failed_outer_ite) { amrex::FArrayBox radEnergySource(indexRange, Physics_Traits::nGroups, amrex::The_Async_Arena()); // cell-centered scalar @@ -1805,7 +1832,8 @@ void QuokkaSimulation::operatorSplitSourceTerms(amrex::Array4::SetRadEnergySource(radEnergySource.array(), indexRange, dx, prob_lo, prob_hi, time + dt); // cell-centered source terms - RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_); + RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_, p_num_failed_coupling, + p_num_failed_dust, p_num_failed_outer_ite); } template diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index cabcf063f..99fe99393 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -205,7 +205,7 @@ template class RadSystem : public HyperbolicSystem::ComputeDustTemperature(double c } } - AMREX_ASSERT_WITH_MESSAGE(ite_td < max_ite_td, "Newton iteration for dust temperature did not converge."); + AMREX_ASSERT_WITH_MESSAGE(ite_td < max_ite_td, "Newton iteration for dust temperature failed to converge."); + if (ite_td >= max_ite_td) { + T_d = -1.0; + } return T_d; } template void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt_radiation, - const int stage, double dustGasCoeff) + const int stage, double dustGasCoeff, int *p_num_failed_coupling, int *p_num_failed_dust, int *p_num_failed_outer_ite) { arrayconst_t &consPrev = consVar; // make read-only array_t &consNew = consVar; @@ -1279,6 +1282,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // cell-centered kernel amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + // make a local reference of p_num_failed + auto p_num_failed_coupling_local = p_num_failed_coupling; + auto p_num_failed_dust_local = p_num_failed_dust; + auto p_num_failed_outer_local = p_num_failed_outer_ite; + const double c = c_light_; const double chat = c_hat_; const double dustGasCoeff_local = dustGasCoeff; @@ -1439,7 +1447,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray F_D{}; const double resid_tol = 1.0e-11; // 1.0e-15; - const int maxIter = 400; + const int maxIter = 50; int n = 0; for (; n < maxIter; ++n) { T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); @@ -1456,8 +1464,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const auto Lambda_gd = sum(Rvec) / (dt * chat / c); T_d = T_gas - Lambda_gd / (dustGasCoeff_local * num_den * num_den * std::sqrt(T_gas)); } + AMREX_ASSERT_WITH_MESSAGE(T_d >= 0., "Dust temperature is negative!"); + if (T_d < 0.0) { + amrex::Gpu::Atomic::Add(p_num_failed_dust_local, 1); + } } - AMREX_ASSERT(T_d >= 0.); fourPiBoverC = ComputeThermalRadiation(T_d, radBoundaries_g_copy); @@ -1677,10 +1688,14 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // } } // END NEWTON-RAPHSON LOOP - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); + AMREX_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); + if (n >= maxIter) { + amrex::Gpu::Atomic::Add(p_num_failed_coupling_local, 1); + } + // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; - AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); - AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); + AMREX_ASSERT(Egas_guess > 0.0); + AMREX_ASSERT(min(EradVec_guess) >= 0.0); if (n > 0) { // calculate kappaF since the temperature has changed @@ -1934,7 +1949,10 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } // end full-step iteration - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ite < max_ite, "AddSourceTerms iteration failed to converge!"); + AMREX_ASSERT_WITH_MESSAGE(ite < max_ite, "AddSourceTerms iteration failed to converge!"); + if (ite >= max_ite) { + amrex::Gpu::Atomic::Add(p_num_failed_outer_local, 1); + } // 4b. Store new radiation energy, gas energy // In the first stage of the IMEX scheme, the hydro quantities are updated by a fraction (defined by