Add 1d version of flux_hllc with normal::AbstractVector (#2276)

* add 1d hllc flux with normal::AbstractVector argument

* add test for 1D hllc flux with normal::AbstractVector

* accomodate non-unit normals

* format

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: jlchan

src/equations/compressible_euler_1d.jl

@@ -816,6 +816,26 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     return SVector(f1, f2, f3)
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that
+# the normal component is incorporated into the numerical flux. 
+#
+# The HLLC flux along a 1D "normal" can be evaluated by scaling the velocity/momentum by 
+# the normal for the 1D HLLC flux, then scaling the resulting momentum flux again. 
+# Moreover, the 2D HLLC flux reduces to this if the normal vector is [n, 0].
+function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector,
+                   equations::CompressibleEulerEquations1D)
+    norm_ = abs(normal_direction[1])
+    normal_direction_unit = normal_direction[1] * inv(norm_)
+
+    # scale the momentum by the normal direction
+    f = flux_hllc(SVector(u_ll[1], normal_direction_unit * u_ll[2], u_ll[3]),
+                  SVector(u_rr[1], normal_direction_unit * u_rr[2], u_rr[3]), 1,
+                  equations)
+
+    # rescale the momentum flux by the normal direction and normalize
+    return SVector(f[1], normal_direction_unit * f[2], f[3]) * norm_
+end
+
 """
     min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)
 

test/test_unit.jl

@@ -1224,6 +1224,25 @@ end
               flux(u, normal_direction, equations)
     end
 
+    # check consistency between 1D and 2D HLLC fluxes
+    u_1d = SVector(1.1, -0.5, 5.5)
+    u_2d = SVector(u_1d[1], u_1d[2], 0.0, u_1d[3])
+    normal_1d = SVector(-0.3)
+    normal_2d = SVector(normal_1d[1], 0.0)
+    equations_1d = CompressibleEulerEquations1D(1.4)
+    equations_2d = CompressibleEulerEquations2D(1.4)
+    flux_1d = flux_hllc(u_1d, u_1d, normal_1d, equations_1d)
+    flux_2d = flux_hllc(u_2d, u_2d, normal_2d, equations_2d)
+    @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d)
+    @test flux_1d ≈ flux_2d[[1, 2, 4]]
+
+    # test when u_ll is not the same as u_rr
+    u_rr_1d = SVector(2.1, 0.3, 0.1)
+    u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3])
+    flux_1d = flux_hllc(u_1d, u_rr_1d, normal_1d, equations_1d)
+    flux_2d = flux_hllc(u_2d, u_rr_2d, normal_2d, equations_2d)
+    @test flux_1d ≈ flux_2d[[1, 2, 4]]
+
     equations = CompressibleEulerEquations3D(1.4)
     u = SVector(1.1, -0.5, 2.34, 2.4, 5.5)