From 9d8aac9c842d05865bef0bb7e50a6922710ed246 Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie.sde@gmail.com>
Date: Fri, 26 Jul 2024 00:38:28 -1000
Subject: [PATCH] Add numerical support for other real types (`traffic_flow`)
 (#2020)

* start

* complete equation

* complete test
---
 src/equations/traffic_flow_lwr_1d.jl | 34 +++++++++++++++-------------
 test/test_type.jl                    | 33 +++++++++++++++++++++++++++
 2 files changed, 51 insertions(+), 16 deletions(-)

diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl
index a4d2613a5c..c41fbb2809 100644
--- a/src/equations/traffic_flow_lwr_1d.jl
+++ b/src/equations/traffic_flow_lwr_1d.jl
@@ -36,11 +36,12 @@ varnames(::typeof(cons2prim), ::TrafficFlowLWREquations1D) = ("car-density",)
 A smooth initial condition used for convergence tests.
 """
 function initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
-    c = 2.0
-    A = 1.0
+    RealT = eltype(x)
+    c = 2
+    A = 1
     L = 1
-    f = 1 / L
-    omega = 2 * pi * f
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
     scalar = c + A * sin(omega * (x[1] - t))
 
     return SVector(scalar)
@@ -55,11 +56,12 @@ Source terms used for convergence tests in combination with
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::TrafficFlowLWREquations1D)
     # Same settings as in `initial_condition`
-    c = 2.0
-    A = 1.0
+    RealT = eltype(x)
+    c = 2
+    A = 1
     L = 1
-    f = 1 / L
-    omega = 2 * pi * f
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
     du = omega * cos(omega * (x[1] - t)) *
          (-1 - equations.v_max * (2 * sin(omega * (x[1] - t)) + 3))
 
@@ -68,21 +70,21 @@ end
 
 # Calculate 1D flux in for a single point
 @inline function flux(u, orientation::Integer, equations::TrafficFlowLWREquations1D)
-    return SVector(equations.v_max * u[1] * (1.0 - u[1]))
+    return SVector(equations.v_max * u[1] * (1 - u[1]))
 end
 
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::TrafficFlowLWREquations1D)
-    λ_max = max(abs(equations.v_max * (1.0 - 2 * u_ll[1])),
-                abs(equations.v_max * (1.0 - 2 * u_rr[1])))
+    λ_max = max(abs(equations.v_max * (1 - 2 * u_ll[1])),
+                abs(equations.v_max * (1 - 2 * u_rr[1])))
 end
 
 # Calculate minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::TrafficFlowLWREquations1D)
-    jac_L = equations.v_max * (1.0 - 2 * u_ll[1])
-    jac_R = equations.v_max * (1.0 - 2 * u_rr[1])
+    jac_L = equations.v_max * (1 - 2 * u_ll[1])
+    jac_R = equations.v_max * (1 - 2 * u_rr[1])
 
     λ_min = min(jac_L, jac_R)
     λ_max = max(jac_L, jac_R)
@@ -96,7 +98,7 @@ end
 end
 
 @inline function max_abs_speeds(u, equations::TrafficFlowLWREquations1D)
-    return (abs(equations.v_max * (1.0 - 2 * u[1])),)
+    return (abs(equations.v_max * (1 - 2 * u[1])),)
 end
 
 # Convert conservative variables to primitive
@@ -106,11 +108,11 @@ end
 @inline cons2entropy(u, equations::TrafficFlowLWREquations1D) = u
 
 # Calculate entropy for a conservative state `cons`
-@inline entropy(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+@inline entropy(u::Real, ::TrafficFlowLWREquations1D) = 0.5f0 * u^2
 @inline entropy(u, equations::TrafficFlowLWREquations1D) = entropy(u[1], equations)
 
 # Calculate total energy for a conservative state `cons`
-@inline energy_total(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+@inline energy_total(u::Real, ::TrafficFlowLWREquations1D) = 0.5f0 * u^2
 @inline energy_total(u, equations::TrafficFlowLWREquations1D) = energy_total(u[1],
                                                                              equations)
 end # @muladd
diff --git a/test/test_type.jl b/test/test_type.jl
index 34fef14fee..6c51460e6d 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -2194,6 +2194,39 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT
         end
     end
+
+    @timed_testset "Traffic Flow LWR 1D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred TrafficFlowLWREquations1D(RealT(1))
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = SVector(one(RealT))
+            orientation = 1
+            c = one(RealT)
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
+                  RealT
+
+            @test eltype(@inferred flux(u, orientation, equations)) == RealT
+
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+            @test typeof(@inferred entropy(c, equations)) == RealT
+            @test typeof(@inferred entropy(u, equations)) == RealT
+            @test typeof(@inferred energy_total(c, equations)) == RealT
+            @test typeof(@inferred energy_total(u, equations)) == RealT
+        end
+    end
 end
 
 end # module