Define K_S for isothermal gamma (#573)

### Description
Currently, `K_S` is undefined in the `computeFlatenningCoefficients`
function for the case of isothermal gamma, because the function to
compute sound speed within `EOS.hpp` that it calls only works if gamma
!= 1.

### Related issues
Fixes #553 for now. In the long run, this wont be needed because we will
move all EOS calculations (including for gamma=1) to EOS.hpp. But for
now, we should apply this patch so that the gamma=1 case is correctly
handled while computing flattening coefficients.

### 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/hydro_system.hpp

@@ -600,8 +600,12 @@ void HydroSystem<problem_t>::ComputeFlatteningCoefficients(amrex::MultiFab const
 
 		// Z is a measure of shock strength (Eq. 76 of Miller & Colella 2002)
 		amrex::GpuArray<Real, nmscalars_> massScalars = RadSystem<problem_t>::ComputeMassScalars(primVar, i, j, k);
-		const double K_S = std::pow(quokka::EOS<problem_t>::ComputeSoundSpeed(primVar(i, j, k, primDensity_index), P, massScalars), 2) *
-				   primVar(i, j, k, primDensity_index);
+		double K_S = std::pow(quokka::EOS<problem_t>::ComputeSoundSpeed(primVar(i, j, k, primDensity_index), P, massScalars), 2) *
+			     primVar(i, j, k, primDensity_index);
+		if constexpr (is_eos_isothermal()) {
+			K_S = primVar(i, j, k, primDensity_index) * cs_iso_ * cs_iso_;
+		}
+
 		const double Z = std::abs(Pplus1 - Pminus1) / K_S;
 
 		// check for converging flow along the normal direction DIR (Eq. 77)