retry hydro update if cooling solve fails (#616)

### Description
This triggers the retry mechanism for the hydro update if the cooling
solve fails. This mirrors the behavior for the chemistry solve (#615).

### Related issues
Closes #372.

### 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.

src/GrackleLikeCooling.hpp

@@ -255,7 +255,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto user_rhs(Real /*t*/, quokka::valar
 	return 0; // success
 }
 
-template <typename problem_t> void computeCooling(amrex::MultiFab &mf, const Real dt_in, grackle_tables &cloudyTables, const Real T_floor)
+template <typename problem_t> auto computeCooling(amrex::MultiFab &mf, const Real dt_in, grackle_tables &cloudyTables, const Real T_floor) -> bool
 {
 	BL_PROFILE("computeCooling()")
 
@@ -309,14 +309,16 @@ template <typename problem_t> void computeCooling(amrex::MultiFab &mf, const Rea
 		});
 	}
 
-	int nmin = nsubstepsMF.min(0);
 	int nmax = nsubstepsMF.max(0);
 	Real navg = static_cast<Real>(nsubstepsMF.sum(0)) / static_cast<Real>(nsubstepsMF.boxArray().numPts());
-	amrex::Print() << fmt::format("\tcooling substeps (per cell): min {}, avg {}, max {}\n", nmin, navg, nmax);
+	amrex::Print() << fmt::format("\tcooling substeps (per cell): avg {}, max {}\n", navg, nmax);
 
+	// check if integration succeeded
 	if (nmax >= maxStepsODEIntegrate) {
-		amrex::Abort("Max steps exceeded in cooling solve!");
+		amrex::Print() << "\t[GrackleLikeCooling] Reaction ODE failure! Retrying hydro update...\n";
+		return false;
 	}
+	return true; // success
 }
 
 void readGrackleData(std::string &grackle_hdf5_file, grackle_tables &cloudyTables);

src/RadhydroSimulation.hpp

@@ -509,21 +509,21 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::addStrangSplit
 template <typename problem_t>
 auto RadhydroSimulation<problem_t>::addStrangSplitSourcesWithBuiltin(amrex::MultiFab &state, int lev, amrex::Real time, amrex::Real dt) -> bool
 {
-
-	// start by assuming chemistry burn is successful.
-	bool burn_success = true; // NOLINT
-
+	// start by assuming cooling integrator is successful.
+	bool cool_success = true;
 	if (enableCooling_ == 1) {
 		// compute cooling
 		if (coolingTableType_ == "grackle") {
-			quokka::GrackleLikeCooling::computeCooling<problem_t>(state, dt, grackleTables_, tempFloor_);
+			cool_success = quokka::GrackleLikeCooling::computeCooling<problem_t>(state, dt, grackleTables_, tempFloor_);
 		} else if (coolingTableType_ == "cloudy_cooling_tools") {
-			quokka::TabulatedCooling::computeCooling<problem_t>(state, dt, cloudyTables_, tempFloor_);
+			cool_success = quokka::TabulatedCooling::computeCooling<problem_t>(state, dt, cloudyTables_, tempFloor_);
 		} else {
 			amrex::Abort("Invalid cooling table type!");
 		}
 	}
 
+	// start by assuming chemistry burn is successful.
+	bool burn_success = true; // NOLINT
 #ifdef PRIMORDIAL_CHEM
 	if (enableChemistry_ == 1) {
 		// compute chemistry
@@ -534,7 +534,7 @@ auto RadhydroSimulation<problem_t>::addStrangSplitSourcesWithBuiltin(amrex::Mult
 	// compute user-specified sources
 	addStrangSplitSources(state, lev, time, dt);
 
-	return burn_success;
+	return (burn_success && cool_success);
 }
 
 template <typename problem_t>
@@ -1036,7 +1036,7 @@ auto RadhydroSimulation<problem_t>::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
+	// check if reactions failed for source terms. If it failed, return false.
 	if (!burn_success_first) {
 		return burn_success_first;
 	}
@@ -1303,7 +1303,7 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 	// do Strang split source terms (second half-step)
 	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
+	// check if we have violated the CFL timestep or reactions failed for source terms
 	return (!isCflViolated(lev, time, dt_lev) && burn_success_second);
 }
 

src/TabulatedCooling.hpp

@@ -310,7 +310,7 @@ template <typename problem_t> auto computeCooling(amrex::MultiFab &mf, const Rea
 
 	// check if integration succeeded
 	if (nmax >= maxStepsODEIntegrate) {
-		amrex::Print() << "\t[CloudyCooling] Reaction ODE failure. Max steps exceeded in cooling solve!\n";
+		amrex::Print() << "\t[CloudyCooling] Reaction ODE failure! Retrying hydro update...\n";
 		return false;
 	}
 	return true; // success