diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 6c6b9dd3721..743d2df870a 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -34,9 +34,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation1D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -49,13 +50,14 @@ in non-periodic domains). function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation1D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -161,7 +163,7 @@ function flux_engquist_osher(u_ll, u_rr, orientation::Int, u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5 * (flux(u_L, orientation, equation) + + return SVector(0.5f0 * (flux(u_L, orientation, equation) + flux(u_R, orientation, equation) - abs(equation.advection_velocity[orientation]) * (u_R - u_L))) end @@ -200,14 +202,16 @@ end @inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer, equations::LinearScalarAdvectionEquation1D) + RealT = eltype(u) a = equations.advection_velocity[1] - return a > 0 ? flux(u, orientation, equations) : zero(u) + return a > 0 ? flux(u, orientation, equations) : SVector(zero(RealT)) end @inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer, equations::LinearScalarAdvectionEquation1D) + RealT = eltype(u) a = equations.advection_velocity[1] - return a < 0 ? flux(u, orientation, equations) : zero(u) + return a < 0 ? flux(u, orientation, equations) : SVector(zero(RealT)) end # Convert conservative variables to primitive @@ -217,11 +221,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation1D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation1D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation1D) energy_total(u[1], equation) end diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl index d90bf0c8793..5e4f8463f52 100644 --- a/src/equations/linear_scalar_advection_2d.jl +++ b/src/equations/linear_scalar_advection_2d.jl @@ -34,8 +34,8 @@ varnames(::typeof(cons2prim), ::LinearScalarAdvectionEquation2D) = ("scalar",) function x_trans_periodic_2d(x, domain_length = SVector(10, 10), center = SVector(0, 0)) x_normalized = x .- center x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5 * domain_length) - - (x_shifted .> 0.5 * domain_length)) .* domain_length + x_offset = ((x_shifted .< -0.5f0 * domain_length) - + (x_shifted .> 0.5f0 * domain_length)) .* domain_length return center + x_shifted + x_offset end @@ -47,9 +47,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x_trans_periodic_2d(x - equation.advection_velocity * t) - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -60,13 +61,14 @@ A smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -280,11 +282,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation2D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation2D) energy_total(u[1], equation) end diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 7b19974eb49..088f934cc3e 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -38,9 +38,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -51,13 +52,14 @@ A smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -82,10 +84,10 @@ A sine wave in the conserved variable. """ function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - scalar = sin(2 * pi * x_trans[1]) * sin(2 * pi * x_trans[2]) * - sin(2 * pi * x_trans[3]) + scalar = sinpi(2 * x_trans[1]) * sinpi(2 * x_trans[2]) * sinpi(2 * x_trans[3]) return SVector(scalar) end @@ -199,11 +201,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation3D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation3D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation3D) energy_total(u[1], equation) end diff --git a/test/test_type.jl b/test/test_type.jl index 7bea104bb86..933d364185e 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -685,6 +685,249 @@ isdir(outdir) && rm(outdir, recursive = true) @test typeof(@inferred density_pressure(u, equations)) == RealT end end + + @timed_testset "Linear Scalar Advection 1D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT)) + orientation = 1 + directions = [1, 2] + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_linear_x(x, t, equations)) == + RealT + + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT + end + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.flux_engquist_osher(u_ll, u_rr, orientation, + equations)) == RealT + + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, + equations))) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Linear Scalar Advection 2D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation2D(RealT(1), RealT(1)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred Trixi.x_trans_periodic_2d(x)) == RealT + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin_sin(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_x_y(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_x(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_y(x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end + end + 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 max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + + for orientation in orientations + @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 + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Linear Scalar Advection 3D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation3D(RealT(1), RealT(1), + RealT(1)) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(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_lax_friedrichs + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_linear_z(x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end + end + 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 max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @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 + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end end end # module