fix numerical flux treatment in subcell scheme

src/GRHD/types.jl

@@ -28,21 +28,24 @@ formulation(::Project{F}) where F = F
 
 struct SubcellDG{T_Mesh<:Mesh{<:SpectralElement}}
   mesh::T_Mesh
-  # ws::FastLapackInterface.LUWs
   Mcube::Array{Float64,3}
   Scube::Array{Float64,3}
+  Bcube::Array{Float64,3}
   M::Matrix{Float64}
   S::Matrix{Float64}
+  B::Matrix{Float64}
+  nflx::Vector{Float64}
 end
 function SubcellDG(mesh::Mesh1d{<:SpectralElement})
   @assert mesh.element.quadr ∈ (:LGL, :LGL_lumped)
   M = similar(mesh.element.invM)
   S = similar(mesh.element.invM)
+  B = similar(mesh.element.invM)
   Mcube = compute_mass_matrix_cube(mesh)
   Scube = compute_stiffness_matrix_cube(mesh)
-  # ws = FastLapackInterface.Workspace(LinearAlgebra.LAPACK.getrf!, M)
-  # return SubcellDG(mesh, ws, Mcube, Scube, M, S)
-  return SubcellDG(mesh, Mcube, Scube, M, S)
+  Bcube = compute_boundary_matrix_cube(mesh)
+  nflx = zeros(Float64, size(B,2))
+  return SubcellDG(mesh, Mcube, Scube, Bcube, M, S, B, nflx)
 end
 
 
@@ -71,20 +74,28 @@ end
 function compute_stiffness_matrix_cube(mesh::Mesh1d{<:SpectralElement})
   @unpack Npts, w, z = mesh.element
   Scube = zeros(Float64, (Npts,Npts,Npts))
+  for i in 1:Npts, j in 1:Npts, n in 1:Npts
+    zn, wi = z[n], w[i]
+    ∂lj_n = dg1d.derivative_lagrange_poly(j, zn, z)
+    Scube[i,j,n] = wi*∂lj_n * (i==n)
+  end
+  return Scube
+end
+function compute_boundary_matrix_cube(mesh::Mesh1d{<:SpectralElement})
+  @unpack Npts, w, z = mesh.element
+  Bcube = zeros(Float64, (Npts,Npts,Npts))
   subz = zeros(Float64, Npts+1)
   subz[1], subz[end] = z[1], z[end]
   for n in 1:Npts-1
     subz[n+1] = z[n]+(z[n+1]-z[n])/2
   end
   for i in 1:Npts, j in 1:Npts, n in 1:Npts
-    zn, wi = z[n], w[i]
     subz_m, subz_p = subz[n], subz[n+1]
     li_m, lj_m = dg1d.lagrange_poly(i, subz_m, z), dg1d.lagrange_poly(j, subz_m, z)
     li_p, lj_p = dg1d.lagrange_poly(i, subz_p, z), dg1d.lagrange_poly(j, subz_p, z)
-    ∂lj_n = dg1d.derivative_lagrange_poly(j, zn, z)
-    Scube[i,j,n] = wi*∂lj_n * (i==n) - li_p*lj_p + li_m*lj_m
+    Bcube[i,j,n] = li_p*lj_p - li_m*lj_m
   end
-  return Scube
+  return Bcube
 end
 
 
@@ -113,10 +124,21 @@ end
     S[i,j] = sij
   end
 end
+@inline function compute_boundary_matrix!(sdg::SubcellDG, vmask)
+  @unpack B, Bcube = sdg
+  @unpack Npts = sdg.mesh.element
+  @turbo for j in 1:Npts, i in 1:Npts
+    bij = 0.0
+    for k in 1:Npts
+      bij += Bcube[i,j,k] * (1.0 - vmask[k])
+    end
+    B[i,j] = bij
+  end
+end
 
 
 function dg1d.compute_rhs_weak_form!(rhs, f, nf, sdg::SubcellDG, mask)
-  @unpack M, S, mesh = sdg
+  @unpack M, S, B, nflx, mesh = sdg
   @unpack w = mesh.element
   @unpack invdetJ = get_static_variables(mesh.cache)
   Npts, K = layout(mesh)
@@ -135,20 +157,18 @@ function dg1d.compute_rhs_weak_form!(rhs, f, nf, sdg::SubcellDG, mask)
     end
     # TODO With lumping approx we can skip mass matrix computation and also divide out
     # the quadrature weight in stiffness matrix computation
-    compute_mass_matrix!(sdg, v_mask)
+    # compute_mass_matrix!(sdg, v_mask)
     compute_stiffness_matrix!(sdg, v_mask)
+    compute_boundary_matrix!(sdg, v_mask)
     v_f    = view(f, cidxs)
     v_invJ = view(invdetJ, cidxs)
-    # mul!(v_rhs, transpose(S), v_f)
-    v_rhs .= transpose(S)*v_f
-    v_rhs[1]   -= nf_lhs[k]
-    v_rhs[end] -= nf_rhs[k]
-    # v_rhs .= M \ v_rhs
-    # v_rhs .= v_rhs ./ diag(M)
-    # v_rhs .*= v_invJ .* v_mask
-    for i in 1:Npts
-      # v_rhs[i] *= v_mask[i] * v_invJ[i] / M[i,i]
-      v_rhs[i] *= (1-v_mask[i]) * v_invJ[i] / w[i]
+    mul!(v_rhs, transpose(S), v_f)
+    nflx[1]   = -nf_lhs[k]
+    @views nflx[2:end-1] .= v_f[2:end-1]
+    nflx[end] = nf_rhs[k]
+    mul!(v_rhs, B, nflx, -1.0, 1.0)
+    @turbo for i in 1:Npts
+      v_rhs[i] *= (1.0 - v_mask[i]) * v_invJ[i] / w[i]
     end
   end
   return