Fix amrex::abort not aborting unconditionally on GPUs (for primordial…

… chem and popiii problems) (#575)

### Description
As described in #543 , `amrex::Abort()` does not unconditionally abort
on GPUs. Following the workaround proposed in #543, I have implemented a
way to circumvent this issue for problems using microphysics' `burn`
(PrimordialChem and PopIIII).

### Related issues
#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.
- [ ] 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>
Co-authored-by: Piyush Sharda <[email protected]>

src/Chemistry.hpp

@@ -31,15 +31,19 @@ AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, Real dt);
 template <typename problem_t> 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<int> 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);
 
-		if (dt < 0) {
-			amrex::Abort("Cannot do chemistry with dt < 0!");
-		}
-
 		amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
 			const Real rho = state(i, j, k, HydroSystem<problem_t>::density_index);
 			const Real xmom = state(i, j, k, HydroSystem<problem_t>::x1Momentum_index);
@@ -60,6 +64,8 @@ template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const R
 			// 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];
@@ -94,7 +100,11 @@ template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const R
 			}
 
 			if (!chemstate.success) {
-				amrex::Abort("VODE integration was unsuccessful!");
+				burn_failed = 1;
+			}
+
+			if (burn_failed) {
+				amrex::Gpu::Atomic::Add(p_num_failed, burn_failed);
 			}
 
 			// ensure positivity and normalize
@@ -142,6 +152,19 @@ template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const R
 				state(i, j, k, HydroSystem<problem_t>::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.
+#endif
+	}
+
+	num_failed = *(d_num_failed.copyToHost());
+
+	burn_success = !num_failed;
+	amrex::ParallelDescriptor::ReduceIntMin(burn_success);
+
+	if (!burn_success) {
+		amrex::Abort("Burn failed in VODE. Aborting.");
 	}
 }