From c485964667aa0ef8324671ba856286cfc8f63287 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 9 Jun 2024 20:54:03 -0700 Subject: [PATCH 1/6] start --- src/equations/linear_scalar_advection_1d.jl | 8 ++++---- src/equations/linear_scalar_advection_2d.jl | 10 +++++----- src/equations/linear_scalar_advection_3d.jl | 6 +++--- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 6c6b9dd3721..2ef1de7dce0 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -52,7 +52,7 @@ function initial_condition_convergence_test(x, t, x_trans = x - equation.advection_velocity * t c = 1.0 - A = 0.5 + A = 0.5f0 L = 2 f = 1 / L omega = 2 * pi * f @@ -161,7 +161,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 @@ -217,11 +217,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..9f027bdada0 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 @@ -63,7 +63,7 @@ function initial_condition_convergence_test(x, t, x_trans = x - equation.advection_velocity * t c = 1.0 - A = 0.5 + A = 0.5f0 L = 2 f = 1 / L omega = 2 * pi * f @@ -280,11 +280,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..ce6d0e6231a 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -54,7 +54,7 @@ function initial_condition_convergence_test(x, t, x_trans = x - equation.advection_velocity * t c = 1.0 - A = 0.5 + A = 0.5f0 L = 2 f = 1 / L omega = 2 * pi * f @@ -199,11 +199,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 From 59070fe72a8d24201ccfd733b1a738e2129a509e Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 10 Jun 2024 22:51:39 -0700 Subject: [PATCH 2/6] complete equations --- src/equations/linear_scalar_advection_1d.jl | 10 ++++++---- src/equations/linear_scalar_advection_2d.jl | 10 ++++++---- src/equations/linear_scalar_advection_3d.jl | 16 ++++++++++------ 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 2ef1de7dce0..6696426c324 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 + 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 diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl index 9f027bdada0..5e4f8463f52 100644 --- a/src/equations/linear_scalar_advection_2d.jl +++ b/src/equations/linear_scalar_advection_2d.jl @@ -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 + 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 diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index ce6d0e6231a..036582fcf72 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 + 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,12 @@ 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 = sin(2 * convert(RealT, pi) * x_trans[1]) * + sin(2 * convert(RealT, pi) * x_trans[2]) * + sin(2 * convert(RealT, pi) * x_trans[3]) return SVector(scalar) end From 22ab4dc88a0a20b6f4ebe5416887b2b48e9ce98a Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 11 Jun 2024 20:32:44 -0700 Subject: [PATCH 3/6] fix zeros --- src/equations/linear_scalar_advection_1d.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 6696426c324..b0879eb504f 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -202,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) : 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) : zero(RealT) end # Convert conservative variables to primitive From 5518bd03c2c4bd1967ab365e1be5b47184c2784a Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Tue, 11 Jun 2024 20:56:30 -0700 Subject: [PATCH 4/6] apply suggestions Co-authored-by: Hendrik Ranocha --- src/equations/linear_scalar_advection_3d.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 036582fcf72..088f934cc3e 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -87,9 +87,7 @@ function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation3D) RealT = eltype(x) x_trans = x - equation.advection_velocity * t - scalar = sin(2 * convert(RealT, pi) * x_trans[1]) * - sin(2 * convert(RealT, pi) * x_trans[2]) * - sin(2 * convert(RealT, pi) * x_trans[3]) + scalar = sinpi(2 * x_trans[1]) * sinpi(2 * x_trans[2]) * sinpi(2 * x_trans[3]) return SVector(scalar) end From 3fca8860f8503ba5f8d8666f36e0bb9c6a303937 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 13 Jun 2024 19:24:19 -1000 Subject: [PATCH 5/6] minor fixes --- src/equations/linear_scalar_advection_1d.jl | 4 ++-- src/equations/linear_scalar_advection_3d.jl | 4 +--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index b0879eb504f..743d2df870a 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -204,14 +204,14 @@ end equations::LinearScalarAdvectionEquation1D) RealT = eltype(u) a = equations.advection_velocity[1] - return a > 0 ? flux(u, orientation, equations) : zero(RealT) + 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(RealT) + return a < 0 ? flux(u, orientation, equations) : SVector(zero(RealT)) end # Convert conservative variables to primitive diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 036582fcf72..088f934cc3e 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -87,9 +87,7 @@ function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation3D) RealT = eltype(x) x_trans = x - equation.advection_velocity * t - scalar = sin(2 * convert(RealT, pi) * x_trans[1]) * - sin(2 * convert(RealT, pi) * x_trans[2]) * - sin(2 * convert(RealT, pi) * x_trans[3]) + scalar = sinpi(2 * x_trans[1]) * sinpi(2 * x_trans[2]) * sinpi(2 * x_trans[3]) return SVector(scalar) end From 23c005034ba2ba58faee3b0ecc32182c67fdc706 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 15 Jun 2024 15:40:29 -1000 Subject: [PATCH 6/6] complete unit tests --- test/test_type.jl | 243 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 243 insertions(+) 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