add LLF fallback solver (#380)

* add LLF fallback solver

* fix typo

src/LLF.hpp

@@ -0,0 +1,48 @@
+#ifndef LLF_HPP_ // NOLINT
+#define LLF_HPP_
+
+#include "AMReX_Extension.H"
+#include "AMReX_GpuQualifiers.H"
+#include <AMReX.H>
+#include <AMReX_REAL.H>
+
+#include "ArrayView.hpp"
+#include "EOS.hpp"
+#include "HydroState.hpp"
+#include "valarray.hpp"
+
+namespace quokka::Riemann
+{
+// Local Lax-Friedrichs (LLF) / Rusanov solver
+template <typename problem_t, int N_scalars, int N_mscalars, int fluxdim>
+AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto LLF(quokka::HydroState<N_scalars, N_mscalars> const &sL, quokka::HydroState<N_scalars, N_mscalars> const &sR)
+    -> quokka::valarray<double, fluxdim>
+{
+	// Toro (Eq. 10.56)
+	const amrex::Real Sp = std::max(std::abs(sL.u) + sL.cs, std::abs(sR.u) + sR.cs);
+
+	/// compute fluxes
+	quokka::valarray<double, fluxdim> U_L = {sL.rho, sL.rho * sL.u, sL.rho * sL.v, sL.rho * sL.w, sL.E, sL.Eint};
+	quokka::valarray<double, fluxdim> U_R = {sR.rho, sR.rho * sR.u, sR.rho * sR.v, sR.rho * sR.w, sR.E, sR.Eint};
+
+	// The remaining components are passive scalars, so just copy them from
+	// x1LeftState and x1RightState into the (left, right) state vectors U_L and
+	// U_R
+	for (int n = 0; n < N_scalars; ++n) {
+		const int nstart = fluxdim - N_scalars;
+		U_L[nstart + n] = sL.scalar[n];
+		U_R[nstart + n] = sR.scalar[n];
+	}
+
+	const quokka::valarray<double, fluxdim> D_L = {0., 1., 0., 0., sL.u, 0.};
+	const quokka::valarray<double, fluxdim> D_R = {0., 1., 0., 0., sR.u, 0.};
+	const quokka::valarray<double, fluxdim> F_L = sL.u * U_L + sL.P * D_L;
+	const quokka::valarray<double, fluxdim> F_R = sR.u * U_R + sR.P * D_R;
+
+	// open the Riemann fan
+	quokka::valarray<double, fluxdim> F = 0.5 * (F_L + F_R) - 0.5 * Sp * (U_R - U_L);
+	return F;
+}
+} // namespace quokka::Riemann
+
+#endif // LLF_HPP_

src/RadhydroSimulation.hpp

@@ -1351,7 +1351,7 @@ void RadhydroSimulation<problem_t>::hydroFluxFunction(amrex::MultiFab const &pri
 	HydroSystem<problem_t>::template FlattenShocks<DIR>(primVar, x1Flat, x2Flat, x3Flat, leftState, rightState, ng_reconstruct, nvars);
 
 	// interface-centered kernel
-	HydroSystem<problem_t>::template ComputeFluxes<DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_);
+	HydroSystem<problem_t>::template ComputeFluxes<RiemannSolver::HLLC, DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_);
 }
 
 template <typename problem_t>
@@ -1401,9 +1401,8 @@ void RadhydroSimulation<problem_t>::hydroFOFluxFunction(amrex::MultiFab const &p
 {
 	// donor-cell reconstruction
 	HydroSystem<problem_t>::template ReconstructStatesConstant<DIR>(primVar, leftState, rightState, ng_reconstruct, nvars);
-
-	// interface-centered kernel
-	HydroSystem<problem_t>::template ComputeFluxes<DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_);
+	// LLF solver
+	HydroSystem<problem_t>::template ComputeFluxes<RiemannSolver::LLF, DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_);
 }
 
 template <typename problem_t> void RadhydroSimulation<problem_t>::swapRadiationState(amrex::MultiFab &stateOld, amrex::MultiFab const &stateNew)

src/hydro_system.hpp

@@ -25,6 +25,7 @@
 #include "ArrayView.hpp"
 #include "EOS.hpp"
 #include "HLLC.hpp"
+#include "LLF.hpp"
 #include "hyperbolic_system.hpp"
 #include "radiation_system.hpp"
 #include "valarray.hpp"
@@ -36,6 +37,10 @@ template <typename problem_t> struct HydroSystem_Traits {
 	static constexpr bool reconstruct_eint = true;
 };
 
+enum class RiemannSolver {
+	HLLC, LLF
+};
+
 /// Class for the Euler equations of inviscid hydrodynamics
 ///
 template <typename problem_t> class HydroSystem : public HyperbolicSystem<problem_t>
@@ -108,7 +113,7 @@ template <typename problem_t> class HydroSystem : public HyperbolicSystem<proble
 
 	static void SyncDualEnergy(amrex::MultiFab &consVar_mf);
 
-	template <FluxDir DIR>
+	template <RiemannSolver RIEMANN, FluxDir DIR>
 	static void ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::MultiFab &x1FaceVel_mf, amrex::MultiFab const &x1LeftState_mf,
 				  amrex::MultiFab const &x1RightState_mf, amrex::MultiFab const &primVar_mf, amrex::Real K_visc);
 
@@ -853,7 +858,7 @@ template <typename problem_t> void HydroSystem<problem_t>::SyncDualEnergy(amrex:
 }
 
 template <typename problem_t>
-template <FluxDir DIR>
+template <RiemannSolver RIEMANN, FluxDir DIR>
 void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::MultiFab &x1FaceVel_mf, amrex::MultiFab const &x1LeftState_mf,
 					   amrex::MultiFab const &x1RightState_mf, amrex::MultiFab const &primVar_mf, const amrex::Real K_visc)
 {
@@ -1029,8 +1034,15 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 		dw = std::min(std::min(dwl, dwr), dw);
 #endif
 
-		// solve the Riemann problem in canonical form
-		quokka::valarray<double, nvar_> F_canonical = quokka::Riemann::HLLC<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR, gamma_, du, dw);
+		// solve the Riemann problem in canonical form (i.e., where the x-dir is the normal direction)
+		quokka::valarray<double, nvar_> F_canonical{};
+
+		if constexpr (RIEMANN == RiemannSolver::HLLC) {
+			F_canonical = quokka::Riemann::HLLC<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR, gamma_, du, dw);
+		} else if constexpr (RIEMANN == RiemannSolver::LLF) {
+			F_canonical = quokka::Riemann::LLF<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR);
+		}
+
 		quokka::valarray<double, nvar_> F = F_canonical;
 
 		// add artificial viscosity