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

…-framework#1971) * start * complete equations * fix zeros * apply suggestions Co-authored-by: Hendrik Ranocha <[email protected]> * minor fixes * complete unit tests --------- Co-authored-by: Hendrik Ranocha <[email protected]>
warisa-r · Jun 16, 2024 · 5398b22 · 5398b22
1 parent 6bb0be5
commit 5398b22
Show file tree

Hide file tree

Showing 4 changed files with 279 additions and 28 deletions.
diff --git 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
@@ -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
@@ -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