fixup dg_rhs in special case when one flux component is nonzero

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

test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml

 [EulerEq]
 id = "kelvin_helmholtz"
+bc = "periodic"
 
 [EquationOfState]
 eos = "idealgas"
@@ -8,12 +9,13 @@ idealgas_gamma = 1.4
 [Mesh]
 range = [ 0.0,1.0, 0.0,1.0 ]
 cfl = 0.1
+# n = [4,4]
 n = [1,1]
-k = [256,256]
+k = [64,64]
 basis = "lgl"
+periodic = [true, true]
 
 [Output]
-every_iteration  = 1
 every_dt  = 0.02
 variables = [ "rho", "qx", "qy", "E", "p", "smoothed_mu" ]
 enable1d  = false
@@ -26,3 +28,5 @@ tend = 2.0
 [HRSC]
 method = "av"
 av_method = "entropy"
+entropy_cmax = 0.1
+entropy_ce = 0.1