Type stability of functions from examples tree_2d_dgsem (#2145)

* Start * Add more * Add more * Add more * Add more * Add more * Add more * Add more * Mark * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * Fix * Apply pi * Fix on second review --------- Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Dec 12, 2024 · ab8ecad · ab8ecad
1 parent 1488a74
commit ab8ecad
Show file tree

Hide file tree

Showing 58 changed files with 575 additions and 506 deletions.
diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl b/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl
@@ -26,9 +26,10 @@ A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
 [`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
 """
 function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)
-    v1_prime = 0.0
-    v2_prime = 0.0
-    p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25)
+    RealT = eltype(x)
+    v1_prime = 0
+    v2_prime = 0
+    p_prime = exp(-log(convert(RealT, 2)) * (x[1]^2 + (x[2] - 25)^2) / 25)
 
     prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)
 

diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl b/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl
@@ -3,18 +3,19 @@ using Trixi
 
 # Oscillating Gaussian-shaped source terms
 function source_terms_gauss(u, x, t, equations::AcousticPerturbationEquations2D)
-    r = 0.1
-    A = 1.0
-    f = 2.0
+    RealT = eltype(u)
+    r = convert(RealT, 0.1)
+    A = 1
+    f = 2
 
     # Velocity sources
-    s1 = 0.0
-    s2 = 0.0
+    s1 = 0
+    s2 = 0
     # Pressure source
-    s3 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sin(2 * pi * f * t)
+    s3 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sinpi(2 * f * t)
 
     # Mean sources
-    s4 = s5 = s6 = s7 = 0.0
+    s4 = s5 = s6 = s7 = 0
 
     return SVector(s1, s2, s3, s4, s5, s6, s7)
 end

diff --git a/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl b/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl
@@ -20,16 +20,17 @@ Initial condition for the monopole in a boundary layer setup, used in combinatio
 [`boundary_condition_monopole`](@ref).
 """
 function initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D)
-    m = 0.3 # Mach number
+    RealT = eltype(x)
+    m = convert(RealT, 0.3) # Mach number
 
-    v1_prime = 0.0
-    v2_prime = 0.0
-    p_prime = 0.0
+    v1_prime = 0
+    v2_prime = 0
+    p_prime = 0
 
     v1_mean = x[2] > 1 ? m : m * (2 * x[2] - 2 * x[2]^2 + x[2]^4)
-    v2_mean = 0.0
-    c_mean = 1.0
-    rho_mean = 1.0
+    v2_mean = 0
+    c_mean = 1
+    rho_mean = 1
 
     prim = SVector(v1_prime, v2_prime, p_prime, v1_mean, v2_mean, c_mean, rho_mean)
 
@@ -48,6 +49,7 @@ with [`initial_condition_monopole`](@ref).
 function boundary_condition_monopole(u_inner, orientation, direction, x, t,
                                      surface_flux_function,
                                      equations::AcousticPerturbationEquations2D)
+    RealT = eltype(u_inner)
     if direction != 3
         error("expected direction = 3, got $direction instead")
     end
@@ -56,9 +58,9 @@ function boundary_condition_monopole(u_inner, orientation, direction, x, t,
     # we use a sinusoidal boundary state for the perturbed variables. For the rest of the -y boundary
     # we set the boundary state to the inner state and multiply the perturbed velocity in the
     # y-direction by -1.
-    if -0.05 <= x[1] <= 0.05 # Monopole
-        v1_prime = 0.0
-        v2_prime = p_prime = sin(2 * pi * t)
+    if RealT(-0.05) <= x[1] <= RealT(0.05) # Monopole
+        v1_prime = 0
+        v2_prime = p_prime = sinpi(2 * t)
 
         prim_boundary = SVector(v1_prime, v2_prime, p_prime, u_inner[4], u_inner[5],
                                 u_inner[6], u_inner[7])

diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl b/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl
@@ -28,7 +28,7 @@ function (indicator::IndicatorAlwaysCoarsen)(u::AbstractArray{<:Any, 4},
     alpha = indicator.cache.alpha
     resize!(alpha, nelements(dg, cache))
 
-    alpha .= -1.0
+    fill!(alpha, -1)
 
     return alpha
 end

diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl b/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl
@@ -28,7 +28,7 @@ function (indicator::IndicatorAlwaysRefine)(u::AbstractArray{<:Any, 4},
     alpha = indicator.cache.alpha
     resize!(alpha, nelements(dg, cache))
 
-    alpha .= 1.0
+    fill!(alpha, 1)
 
     return alpha
 end

diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl
@@ -19,16 +19,17 @@ end
 function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4},
                                                    mesh, equations, dg, cache;
                                                    t, kwargs...)
+    RealT = eltype(u)
     mesh = indicator.cache.mesh
     alpha = indicator.cache.alpha
     resize!(alpha, nelements(dg, cache))
 
     #Predict the theoretical center.
-    advection_velocity = (0.2, -0.7)
+    advection_velocity = (convert(RealT, 0.2), convert(RealT, -0.7))
     center = t .* advection_velocity
 
     inner_distance = 1
-    outer_distance = 1.85
+    outer_distance = convert(RealT, 1.85)
 
     #Iterate over all elements
     for element in eachindex(alpha)
@@ -38,16 +39,16 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4},
 
         #The geometric shape of the amr should be preserved when the base_level is increased.
         #This is done by looking at the original coordinates of each cell.
-        cell_coordinates = original_coordinates(coordinates, 5 / 8)
+        cell_coordinates = original_coordinates(coordinates, 0.625f0)
         cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
         if cell_distance < (inner_distance + outer_distance) / 2
-            cell_coordinates = original_coordinates(coordinates, 5 / 16)
+            cell_coordinates = original_coordinates(coordinates, 0.3125f0)
             cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
         end
 
         #Set alpha according to cells position inside the circles.
         target_level = (cell_distance < inner_distance) + (cell_distance < outer_distance)
-        alpha[element] = target_level / 2
+        alpha[element] = convert(RealT, target_level) / 2
     end
     return alpha
 end

diff --git a/examples/tree_2d_dgsem/elixir_advection_callbacks.jl b/examples/tree_2d_dgsem/elixir_advection_callbacks.jl
@@ -130,7 +130,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2cons)
 
-example_callback = TrixiExtensionExample.ExampleStepCallback(message = "안녕하세요?")
+example_callback = TrixiExtensionExample.ExampleStepCallback(message = "Initializing callback")
 
 stepsize_callback = StepsizeCallback(cfl = 1.6)
 

diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl
@@ -25,14 +25,15 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 function initial_condition_diffusive_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
 
     nu = diffusivity()
-    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)) * exp(-2 * nu * omega^2 * t)
     return SVector(scalar)
 end

diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -28,14 +28,15 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 #   to numerical partial differential equations.
 #   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
 function initial_condition_eriksson_johnson(x, t, equations)
+    RealT = eltype(x)
     l = 4
     epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt
     lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
     lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
-    r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
-    s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+    r1 = (1 + sqrt(1 + 4 * convert(RealT, pi)^2 * epsilon^2)) / (2 * epsilon)
+    s1 = (1 - sqrt(1 + 4 * convert(RealT, pi)^2 * epsilon^2)) / (2 * epsilon)
     u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) +
-        cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
+        cospi(x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
     return SVector{1}(u)
 end
 initial_condition = initial_condition_eriksson_johnson

diff --git a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl
@@ -13,18 +13,19 @@ equations = CompressibleEulerEquations2D(gamma)
 #   https://tinyurl.com/c76fjtx4
 # Mach = 2000 jet
 function initial_condition_astro_jet(x, t, equations::CompressibleEulerEquations2D)
+    RealT = eltype(x)
     @unpack gamma = equations
-    rho = 0.5
+    rho = 0.5f0
     v1 = 0
     v2 = 0
-    p = 0.4127
+    p = convert(RealT, 0.4127)
     # add inflow for t>0 at x=-0.5
     # domain size is [-0.5,+0.5]^2
-    if (t > 0) && (x[1] ≈ -0.5) && (abs(x[2]) < 0.05)
-        rho = 5
+    if (t > 0) && (x[1] ≈ -0.5f0) && (abs(x[2]) < RealT(0.05))
+        rho = 5.0f0
         v1 = 800 # about Mach number Ma = 2000
         v2 = 0
-        p = 0.4127
+        p = convert(RealT, 0.4127)
     end
     return prim2cons(SVector(rho, v1, v2, p), equations)
 end

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl
@@ -18,18 +18,19 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0, 0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
     phi = atan(y_norm, x_norm)
     sin_phi, cos_phi = sincos(phi)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
-    p = r > 0.5 ? 1.0E-3 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : RealT(1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, v2, p), equations)
 end

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl
@@ -18,18 +18,19 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0, 0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
     phi = atan(y_norm, x_norm)
     sin_phi, cos_phi = sincos(phi)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
-    p = r > 0.5 ? 1.0E-3 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : RealT(1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, v2, p), equations)
 end

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl
@@ -18,18 +18,19 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0, 0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
     phi = atan(y_norm, x_norm)
     sin_phi, cos_phi = sincos(phi)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
-    p = r > 0.5 ? 1.0E-3 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : RealT(1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, v2, p), equations)
 end

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl
@@ -18,18 +18,19 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0, 0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
     phi = atan(y_norm, x_norm)
     sin_phi, cos_phi = sincos(phi)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
-    p = r > 0.5 ? 1.0E-3 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : RealT(1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, v2, p), equations)
 end

diff --git a/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl b/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl
@@ -23,17 +23,18 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
     # resolution 128^2, 256^2
     # domain size is [-20.0,20.0]^2
     # gamma = 5/3 for this test case
-    R = 1.0 # radius of the blob
+    RealT = eltype(x)
+    R = 1 # radius of the blob
     # background density
-    dens0 = 1.0
-    Chi = 10.0 # density contrast
+    dens0 = 1
+    Chi = convert(RealT, 10) # density contrast
     # reference time of characteristic growth of KH instability equal to 1.0
-    tau_kh = 1.0
-    tau_cr = tau_kh / 1.6 # crushing time
+    tau_kh = 1
+    tau_cr = tau_kh / convert(RealT, 1.6) # crushing time
     # determine background velocity
     velx0 = 2 * R * sqrt(Chi) / tau_cr
-    vely0 = 0.0
-    Ma0 = 2.7 # background flow Mach number Ma=v/c
+    vely0 = 0
+    Ma0 = convert(RealT, 2.7) # background flow Mach number Ma=v/c
     c = velx0 / Ma0 # sound speed
     # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
     p0 = c * c * dens0 / equations.gamma
@@ -45,9 +46,10 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
     slope = 2
     # density blob
     dens = dens0 +
-           (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
+           (Chi - 1) * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
     # velocity blob is zero
-    velx = velx0 - velx0 * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
+    velx = velx0 -
+           velx0 * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
     return prim2cons(SVector(dens, velx, vely0, p0), equations)
 end
 initial_condition = initial_condition_blob

diff --git a/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl b/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl
@@ -23,17 +23,18 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
     # resolution 128^2, 256^2
     # domain size is [-20.0,20.0]^2
     # gamma = 5/3 for this test case
-    R = 1.0 # radius of the blob
+    RealT = eltype(x)
+    R = 1 # radius of the blob
     # background density
-    dens0 = 1.0
-    Chi = 10.0 # density contrast
+    dens0 = 1
+    Chi = convert(RealT, 10) # density contrast
     # reference time of characteristic growth of KH instability equal to 1.0
-    tau_kh = 1.0
-    tau_cr = tau_kh / 1.6 # crushing time
+    tau_kh = 1
+    tau_cr = tau_kh / convert(RealT, 1.6) # crushing time
     # determine background velocity
     velx0 = 2 * R * sqrt(Chi) / tau_cr
-    vely0 = 0.0
-    Ma0 = 2.7 # background flow Mach number Ma=v/c
+    vely0 = 0
+    Ma0 = convert(RealT, 2.7) # background flow Mach number Ma=v/c
     c = velx0 / Ma0 # sound speed
     # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
     p0 = c * c * dens0 / equations.gamma
@@ -45,9 +46,10 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
     slope = 2
     # density blob
     dens = dens0 +
-           (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
+           (Chi - 1) * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
     # velocity blob is zero
-    velx = velx0 - velx0 * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
+    velx = velx0 -
+           velx0 * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1)))
     return prim2cons(SVector(dens, velx, vely0, p0), equations)
 end
 initial_condition = initial_condition_blob