Skip to content

Commit

Permalink
Make assertion work in the exchange step and dust temperature solutio…
Browse files Browse the repository at this point in the history 
…n 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>
  • Loading branch information
@chongchonghe @pre-commit-ci
chongchonghe and pre-commit-ci[bot] authored Aug 19, 2024
1 parent 21ba4c2 commit 8668640
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 14 deletions.
38 changes: 33 additions & 5 deletions src/QuokkaSimulation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,8 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob

void operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, amrex::Real time, double dt, int stage,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi);
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi, int *num_failed_coupling, int *num_failed_dust,
int *p_num_failed_outer_ite);

auto computeRadiationFluxes(amrex::Array4<const amrex::Real> const &consVar, const amrex::Box &indexRange, int nvars,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx)
Expand Down Expand Up @@ -1615,8 +1616,16 @@ void QuokkaSimulation<problem_t>::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<int> num_failed_coupling({0});
amrex::Gpu::Buffer<int> num_failed_dust({0});
amrex::Gpu::Buffer<int> 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);
Expand All @@ -1625,7 +1634,8 @@ void QuokkaSimulation<problem_t>::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);
}
}

Expand All @@ -1642,7 +1652,23 @@ void QuokkaSimulation<problem_t>::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_
Expand Down Expand Up @@ -1794,7 +1820,8 @@ template <typename problem_t>
void QuokkaSimulation<problem_t>::operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, const amrex::Real time,
const double dt, const int stage, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi)
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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<problem_t>::nGroups,
amrex::The_Async_Arena()); // cell-centered scalar
Expand All @@ -1805,7 +1832,8 @@ void QuokkaSimulation<problem_t>::operatorSplitSourceTerms(amrex::Array4<amrex::
RadSystem<problem_t>::SetRadEnergySource(radEnergySource.array(), indexRange, dx, prob_lo, prob_hi, time + dt);

// cell-centered source terms
RadSystem<problem_t>::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_);
RadSystem<problem_t>::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_, p_num_failed_coupling,
p_num_failed_dust, p_num_failed_outer_ite);
}

template <typename problem_t>
Expand Down
36 changes: 27 additions & 9 deletions src/radiation/radiation_system.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@ template <typename problem_t> class RadSystem : public HyperbolicSystem<problem_
amrex::Real time);

static void AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt, int stage,
double dustGasCoeff);
double dustGasCoeff, int *num_failed_coupling, int *num_failed_dust, int *p_num_failed_outer_ite);

static void balanceMatterRadiation(arrayconst_t &consPrev, array_t &consNew, amrex::Box const &indexRange);

Expand Down Expand Up @@ -1255,13 +1255,16 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::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 <typename problem_t>
void RadSystem<problem_t>::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;
Expand All @@ -1279,6 +1282,11 @@ void RadSystem<problem_t>::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;
Expand Down Expand Up @@ -1439,7 +1447,7 @@ void RadSystem<problem_t>::AddSourceTerms(array_t &consVar, arrayconst_t &radEne
quokka::valarray<double, nGroups_> 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<problem_t>::ComputeTgasFromEint(rho, Egas_guess, massScalars);
Expand All @@ -1456,8 +1464,11 @@ void RadSystem<problem_t>::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);

Expand Down Expand Up @@ -1677,10 +1688,14 @@ void RadSystem<problem_t>::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
Expand Down Expand Up @@ -1934,7 +1949,10 @@ void RadSystem<problem_t>::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
Expand Down

0 comments on commit 8668640

Please sign in to comment.