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/dg_rhs.jl

@@ -73,29 +73,37 @@ function compute_rhs_weak_form!(rhs, fx, fy::Real, nf, mesh::Mesh2d)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   Ml_down, Ml_up = Ml_lhs, Ml_rhs
   Nx, Ny = Npts, Npts
-  @inbounds for k = 1:n_cells(mesh.tree)
+  @inbounds for kk = 1:n_cells(mesh.tree)
 
     # bulk
-    idxs_bulk = cellindices(mesh, k)
+    idxs_bulk = cellindices(mesh, kk)
     v_rhs     = view(rhs,  idxs_bulk)
     v_fx      = view(fx,   idxs_bulk)
     v_dxdX    = view(dxdX, idxs_bulk)
+    v_dydY    = view(dydY, idxs_bulk)
 
     # faces
-    idxs_lhs, idxs_rhs, _, _ = faceindices(mesh, k)
+    idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
 
     v_nf_lhs  = view(nf, idxs_lhs)
     v_nf_rhs  = view(nf, idxs_rhs)
+    v_nf_down = view(nf, idxs_down)
+    v_nf_up   = view(nf, idxs_up)
 
     # TODO Benchmark and see whether LoopVectorization can used here.
     @inbounds for j=1:Ny, i=1:Nx
       rhsij = 0
       dxdXij = v_dxdX[i,j]
+      dydYij = v_dydY[i,j]
       for k = 1:Nx
         rhsij += MDM[i,k] * v_fx[k,j]
       end
       v_rhs[i,j] = rhsij * dxdXij
       v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij
+      # TODO Can we really neglect num fluxes in orthogonal direction in general?
+      #      Or we need to adjust for triangular meshes?
+      # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij +
+      #               (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
     end
 
   end
@@ -107,29 +115,37 @@ function compute_rhs_weak_form!(rhs, fx::Real, fy, nf, mesh::Mesh2d)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   Ml_down, Ml_up = Ml_lhs, Ml_rhs
   Nx, Ny = Npts, Npts
-  @inbounds for k = 1:n_cells(mesh.tree)
+  @inbounds for kk = 1:n_cells(mesh.tree)
 
     # bulk
-    idxs_bulk = cellindices(mesh, k)
+    idxs_bulk = cellindices(mesh, kk)
     v_rhs     = view(rhs,  idxs_bulk)
     v_fy      = view(fy,   idxs_bulk)
+    v_dxdX    = view(dxdX, idxs_bulk)
     v_dydY    = view(dydY, idxs_bulk)
 
     # faces
-    _, _, idxs_down, idxs_up = faceindices(mesh, k)
+    idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
 
+    v_nf_lhs  = view(nf, idxs_lhs)
+    v_nf_rhs  = view(nf, idxs_rhs)
     v_nf_down = view(nf, idxs_down)
     v_nf_up   = view(nf, idxs_up)
 
     # TODO Benchmark and see whether LoopVectorization can used here.
     @inbounds for j=1:Ny, i=1:Nx
       rhsij = 0
+      dxdXij = v_dxdX[i,j]
       dydYij = v_dydY[i,j]
       for k = 1:Ny
         rhsij += MDM[j,k] * v_fy[i,k]
       end
-      v_rhs[i,j] += rhsij * dydYij
+      v_rhs[i,j] = rhsij * dydYij
       v_rhs[i,j] -= (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
+      # TODO Can we really neglect num fluxes in orthogonal direction in general?
+      #      Or we need to adjust for triangular meshes?
+      # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij +
+      #               (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
     end
 
   end