Add numerical support of other real types (linearized_euler) (trixi…

…-framework#2002) * start * complete equations * unit test 1D * unit test 2D * unit test 3D
DanielDoehring · Aug 15, 2024 · 19ab082 · 19ab082
1 parent 0cc13dd
commit 19ab082
Show file tree

Hide file tree

Showing 2 changed files with 167 additions and 13 deletions.
diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl
@@ -204,14 +204,14 @@ The diagonalization of the flux matrix can be found in
         lambda1_p = positive_part(lambda1)
         lambda2_p = positive_part(lambda2)
         lambda3_p = positive_part(lambda3)
-        lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
-        lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+        lambda2p3_half_p = 0.5f0 * (lambda2_p + lambda3_p)
+        lambda3m2_half_p = 0.5f0 * (lambda3_p - lambda2_p)
 
         lambda1_m = negative_part(lambda1)
         lambda2_m = negative_part(lambda2)
         lambda3_m = negative_part(lambda3)
-        lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
-        lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+        lambda2p3_half_m = 0.5f0 * (lambda2_m + lambda3_m)
+        lambda3m2_half_m = 0.5f0 * (lambda3_m - lambda2_m)
 
         f1p = (lambda1_p * rho_prime_ll +
                lambda3m2_half_p / c_mean_global * rho_mean_global * v1_prime_ll +
@@ -244,14 +244,14 @@ The diagonalization of the flux matrix can be found in
         lambda1_p = positive_part(lambda1)
         lambda2_p = positive_part(lambda2)
         lambda3_p = positive_part(lambda3)
-        lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
-        lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+        lambda2p3_half_p = 0.5f0 * (lambda2_p + lambda3_p)
+        lambda3m2_half_p = 0.5f0 * (lambda3_p - lambda2_p)
 
         lambda1_m = negative_part(lambda1)
         lambda2_m = negative_part(lambda2)
         lambda3_m = negative_part(lambda3)
-        lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
-        lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+        lambda2p3_half_m = 0.5f0 * (lambda2_m + lambda3_m)
+        lambda3m2_half_m = 0.5f0 * (lambda3_m - lambda2_m)
 
         f1p = (lambda1_p * rho_prime_ll +
                lambda3m2_half_p / c_mean_global * rho_mean_global * v2_prime_ll +
@@ -304,14 +304,14 @@ end
     lambda1_p = positive_part(lambda1)
     lambda2_p = positive_part(lambda2)
     lambda3_p = positive_part(lambda3)
-    lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
-    lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+    lambda2p3_half_p = 0.5f0 * (lambda2_p + lambda3_p)
+    lambda3m2_half_p = 0.5f0 * (lambda3_p - lambda2_p)
 
     lambda1_m = negative_part(lambda1)
     lambda2_m = negative_part(lambda2)
     lambda3_m = negative_part(lambda3)
-    lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
-    lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+    lambda2p3_half_m = 0.5f0 * (lambda2_m + lambda3_m)
+    lambda3m2_half_m = 0.5f0 * (lambda3_m - lambda2_m)
 
     f1p = (lambda1_p * rho_prime_ll +
            lambda3m2_half_p / c_mean_global * rho_mean_global * v_prime_normal_ll +

diff --git a/test/test_type.jl b/test/test_type.jl
@@ -176,7 +176,7 @@ isdir(outdir) && rm(outdir, recursive = true)
 
     @timed_testset "Compressible Euler 2D" begin
         for RealT in (Float32, Float64)
-            # set gamma = 2 for the coupling convergence test\
+            # set gamma = 2 for the coupling convergence test
             equations = @inferred CompressibleEulerEquations2D(RealT(2))
 
             x = SVector(zero(RealT), zero(RealT))
@@ -1772,6 +1772,160 @@ isdir(outdir) && rm(outdir, recursive = true)
         end
     end
 
+    @timed_testset "Linearized Euler 1D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred LinearizedEulerEquations1D(v_mean_global = RealT(0),
+                                                             c_mean_global = RealT(1),
+                                                             rho_mean_global = RealT(1))
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT))
+            orientation = 1
+            directions = [1, 2]
+
+            surface_flux_function = flux_hll
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+
+            for direction in directions
+                @test eltype(@inferred boundary_condition_wall(u_inner, orientation,
+                                                               direction, x, t,
+                                                               surface_flux_function,
+                                                               equations)) == RealT
+            end
+
+            @test eltype(@inferred flux(u, orientation, equations)) == RealT
+
+            @test typeof(@inferred Trixi.max_abs_speeds(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 cons2prim(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+        end
+    end
+
+    @timed_testset "Linearized Euler 2D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred LinearizedEulerEquations2D(v_mean_global = (RealT(0),
+                                                                              RealT(0)),
+                                                             c_mean_global = RealT(1),
+                                                             rho_mean_global = RealT(1))
+
+            x = SVector(zero(RealT), zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT),
+                                                one(RealT))
+            orientations = [1, 2]
+            directions = [1, 2, 3, 4]
+            normal_direction = SVector(one(RealT), zero(RealT))
+
+            surface_flux_function = flux_hll
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+
+            for orientation in orientations
+                for direction in directions
+                    @test eltype(@inferred boundary_condition_wall(u_inner, orientation,
+                                                                   direction, x, t,
+                                                                   surface_flux_function,
+                                                                   equations)) == RealT
+                end
+
+                @test eltype(@inferred flux(u, orientation, equations)) == RealT
+                @test eltype(@inferred flux_godunov(u_ll, u_rr, 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
+            end
+
+            @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
+            @test eltype(@inferred flux_godunov(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+
+            @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+        end
+    end
+
+    @timed_testset "Linearized Euler 3D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred LinearizedEulerEquations3D(v_mean_global = (RealT(0),
+                                                                              RealT(0),
+                                                                              RealT(0)),
+                                                             c_mean_global = RealT(1),
+                                                             rho_mean_global = RealT(1))
+
+            x = SVector(zero(RealT), zero(RealT), zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT),
+                                                one(RealT), one(RealT))
+            orientations = [1, 2, 3]
+            directions = [1, 2, 3, 4, 5, 6]
+            normal_direction = SVector(one(RealT), zero(RealT), zero(RealT))
+
+            surface_flux_function = flux_hll
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+
+            for orientation in orientations
+                for direction in directions
+                    @test eltype(@inferred boundary_condition_wall(u_inner, orientation,
+                                                                   direction, x, t,
+                                                                   surface_flux_function,
+                                                                   equations)) == RealT
+                end
+
+                @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
+            end
+
+            @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
+
+            @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+        end
+    end
+
     @timed_testset "Shallow Water 1D" begin
         for RealT in (Float32, Float64)
             equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81))