From 1e1f64314d49386ef1a671c0cb5dfda13a048c31 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Thu, 13 Feb 2025 01:49:44 -0600
Subject: [PATCH] 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 <ranocha@users.noreply.github.com>
---
 src/equations/compressible_euler_1d.jl | 20 ++++++++++++++++++++
 test/test_unit.jl                      | 19 +++++++++++++++++++
 2 files changed, 39 insertions(+)

diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
index 542a5b7039f..28572cdf9ee 100644
--- a/src/equations/compressible_euler_1d.jl
+++ b/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)
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index d71a9daea10..2613a448ad9 100644
--- a/test/test_unit.jl
+++ b/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)