Make EulerEq2d run with Kelvin Helmholtz test

• fixed dg weak form computation in 2d for special case when one component of the flux is zero
• fixed/implemented flux and ldg flux equations for 2d
• converted all max calls to absmax calls in equations.jl
• extracted numerical utilities from utils.jl into numutils.jl

The IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml test now runs with polynomial orders n=1 (up to n=4; but not fully tested). Unfortunately, atm this runs pretty slow, so it is not yet added as ref test, but goal is to make it one in the future.

src/EulerEq/rhs.jl

@@ -235,6 +235,7 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
 
   ## solve auxiliary equation: q + ∂x u = 0
   broadcast_faces_2!(ldg_nflux, equation, mesh)
+  # TODO Need to revise Bdry conditions
   # broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
   compute_rhs_weak_form!(ldg_rho_x, rho, 0,   nflx_rho_x, mesh)
   compute_rhs_weak_form!(ldg_rho_y, 0,   rho, nflx_rho_y, mesh)
@@ -248,6 +249,7 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
   ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf, equation, mesh)
+  # TODO Need to revise Bdry conditions
   # broadcast_bdry_2!(bdryllf, equation, mesh)
 
   broadcast_volume_2!(av_flux_2d, equation, mesh)