Merge branch 'main' into example2

examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl

@@ -8,8 +8,8 @@ equations = TrafficFlowLWREquations1D()
 
 solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis))
 
-coordinates_min = (-1.0,) # minimum coordinate
-coordinates_max = (1.0,) # maximum coordinate
+coordinates_min = -1.0 # minimum coordinate
+coordinates_max = 1.0 # maximum coordinate
 cells_per_dimension = (64,)
 
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,

examples/tree_1d_dgsem/elixir_advection_diffusion.jl

@@ -22,10 +22,12 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
                 n_cells_max = 30_000, # set maximum capacity of tree data structure
                 periodicity = true)
 
-function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)),
+                          center = SVector(oftype(x[1], 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)) .*
+    x_offset = ((x_shifted .< -0.5f0 * domain_length) -
+                (x_shifted .> 0.5f0 * domain_length)) .*
                domain_length
     return center + x_shifted + x_offset
 end
@@ -37,11 +39,9 @@ function initial_condition_diffusive_convergence_test(x, t,
     x_trans = x_trans_periodic(x - equation.advection_velocity * t)
 
     nu = diffusivity()
-    c = 0.0
-    A = 1.0
-    L = 2
-    f = 1 / L
-    omega = 1.0
+    c = 0
+    A = 1
+    omega = 1
     scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)
     return SVector(scalar)
 end

examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl

@@ -8,8 +8,9 @@ using Trixi
 equations = InviscidBurgersEquation1D()
 
 function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
+    RealT = eltype(x)
     k = 1
-    2 + sinpi(k * (x[1] - 0.7)) |> SVector
+    SVector(2 + sinpi(k * (x[1] - convert(RealT, 0.7))))
 end
 
 volume_flux = flux_ec

examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl

@@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
 
 # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan.
 function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D)
-    scalar = x[1] < 0.5 ? 0.5 : 1.5
+    RealT = eltype(x)
+    scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0)
 
     return SVector(scalar)
 end

examples/tree_1d_dgsem/elixir_burgers_shock.jl

@@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
 
 # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed.
 function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D)
-    scalar = x[1] < 0.5 ? 1.5 : 0.5
+    RealT = eltype(x)
+    scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0)
 
     return SVector(scalar)
 end

examples/tree_1d_dgsem/elixir_euler_blast_wave.jl

@@ -18,7 +18,8 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
     # The following code is equivalent to
@@ -28,9 +29,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
     cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_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
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, p), equations)
 end

examples/tree_1d_dgsem/elixir_euler_positivity.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

@@ -18,10 +18,11 @@ A discontinuous initial condition taken from
 """
 function initial_condition_discontinuity(x, t,
                                          equations::CompressibleEulerEquationsQuasi1D)
-    rho = (x[1] < 0) ? 3.4718 : 2.0
-    v1 = (x[1] < 0) ? -2.5923 : -3.0
-    p = (x[1] < 0) ? 5.7118 : 2.639
-    a = (x[1] < 0) ? 1.0 : 1.5
+    RealT = eltype(x)
+    rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0)
+    v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0)
+    p = (x[1] < 0) ? RealT(5.7118) : RealT(2.639)
+    a = (x[1] < 0) ? 1.0f0 : 1.5f0
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end

examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl

@@ -15,10 +15,11 @@ equations = CompressibleEulerEquationsQuasi1D(1.4)
 # refinement level is changed the initial condition below may need changed as well to
 # ensure that the discontinuities lie on an element interface.
 function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D)
-    v1 = 0.1
-    rho = 2.0 + 0.1 * x[1]
-    p = 3.0
-    a = 2.0 + x[1]
+    RealT = eltype(x)
+    v1 = convert(RealT, 0.1)
+    rho = 2 + convert(RealT, 0.1) * x[1]
+    p = 3
+    a = 2 + x[1]
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end

examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl

@@ -18,20 +18,21 @@ A multicomponent two interacting blast wave test taken from
 """
 function initial_condition_two_interacting_blast_waves(x, t,
                                                        equations::CompressibleEulerMulticomponentEquations1D)
-    rho1 = 0.5 * x[1]^2
-    rho2 = 0.5 * (sin(20 * x[1]))^2
+    RealT = eltype(x)
+    rho1 = 0.5f0 * x[1]^2
+    rho2 = 0.5f0 * (sin(20 * x[1]))^2
     rho3 = 1 - rho1 - rho2
 
     prim_rho = SVector{3, real(equations)}(rho1, rho2, rho3)
 
-    v1 = 0.0
+    v1 = 0
 
-    if x[1] <= 0.1
-        p = 1000
-    elseif x[1] < 0.9
-        p = 0.01
+    if x[1] <= RealT(0.1)
+        p = convert(RealT, 1000)
+    elseif x[1] < RealT(0.9)
+        p = convert(RealT, 0.01)
     else
-        p = 100
+        p = convert(RealT, 100)
     end
 
     prim_other = SVector{2, real(equations)}(v1, p)

examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl

@@ -19,12 +19,13 @@ A non-priodic harmonic function used in combination with
 function initial_condition_harmonic_nonperiodic(x, t,
                                                 equations::HyperbolicDiffusionEquations1D)
     # elliptic equation: -νΔϕ = f
-    if t == 0.0
-        phi = 5.0
-        q1 = 0.0
+    RealT = eltype(x)
+    if t == 0
+        phi = convert(RealT, 5)
+        q1 = zero(RealT)
     else
         A = 3
-        B = exp(1)
+        B = exp(one(RealT))
         phi = A + B * x[1]
         q1 = B
     end

examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl

@@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 # that is advected to left with v - c and to the right with v + c.
 # Correspondigly, the bump splits in half.
 function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D)
-    v1_prime = 0.0
+    v1_prime = 0
     rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25)
     return SVector(rho_prime, v1_prime, p_prime)
 end

examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl

@@ -20,9 +20,10 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 # Excite the electric field which causes a standing wave.
 # The solution is an undamped exchange between electric and magnetic energy.
 function initial_condition_E_excitation(x, t, equations::MaxwellEquations1D)
+    RealT = eltype(x)
     c = equations.speed_of_light
-    E = -c * sin(2 * pi * x[1])
-    B = 0.0
+    E = -c * sinpi(2 * x[1])
+    B = 0
 
     return SVector(E, B)
 end

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl

@@ -18,14 +18,15 @@ MHD extension of the Sod shock tube. Taken from Section V of the article
 """
 function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
     # domain must be set to [0, 1], γ = 2, final time = 0.12
-    rho = x[1] < 0.5 ? 1.0 : 0.125
-    v1 = 0.0
-    v2 = 0.0
-    v3 = 0.0
-    p = x[1] < 0.5 ? 1.0 : 0.1
-    B1 = 0.75
-    B2 = x[1] < 0.5 ? 1.0 : -1.0
-    B3 = 0.0
+    RealT = eltype(x)
+    rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0
+    v1 = 0
+    v2 = 0
+    v3 = 0
+    p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1)
+    B1 = 0.75f0
+    B2 = x[1] < 0.5f0 ? 1 : -1
+    B3 = 0
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end
 initial_condition = initial_condition_briowu_shock_tube

examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl

@@ -23,14 +23,15 @@ present in the one dimensional MHD equations. It is the second test from Section
 """
 function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
     # domain must be set to [0, 1], γ = 5/3, final time = 0.2
-    rho = x[1] <= 0.5 ? 1.08 : 1.0
-    v1 = x[1] <= 0.5 ? 1.2 : 0.0
-    v2 = x[1] <= 0.5 ? 0.01 : 0.0
-    v3 = x[1] <= 0.5 ? 0.5 : 0.0
-    p = x[1] <= 0.5 ? 0.95 : 1.0
-    inv_sqrt4pi = 1.0 / sqrt(4 * pi)
+    RealT = eltype(x)
+    rho = x[1] <= 0.5f0 ? RealT(1.08) : one(RealT)
+    v1 = x[1] <= 0.5f0 ? RealT(1.2) : zero(RealT)
+    v2 = x[1] <= 0.5f0 ? RealT(0.01) : zero(RealT)
+    v3 = x[1] <= 0.5f0 ? 0.5f0 : 0.0f0
+    p = x[1] <= 0.5f0 ? RealT(0.95) : one(RealT)
+    inv_sqrt4pi = 1 / sqrt(4 * convert(RealT, pi))
     B1 = 2 * inv_sqrt4pi
-    B2 = x[1] <= 0.5 ? 3.6 * inv_sqrt4pi : 4.0 * inv_sqrt4pi
+    B2 = x[1] <= 0.5f0 ? RealT(3.6) * inv_sqrt4pi : 4 * inv_sqrt4pi
     B3 = B1
 
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)

examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl

@@ -19,15 +19,16 @@ Taken from Section 4.1 of
 function initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
     # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
     # initial shock location is taken to be at x = -4
-    x_0 = -4.0
-    rho = x[1] <= x_0 ? 3.5 : 1.0 + 0.2 * sin(5.0 * x[1])
-    v1 = x[1] <= x_0 ? 5.8846 : 0.0
-    v2 = x[1] <= x_0 ? 1.1198 : 0.0
-    v3 = 0.0
-    p = x[1] <= x_0 ? 42.0267 : 1.0
-    B1 = 1.0
-    B2 = x[1] <= x_0 ? 3.6359 : 1.0
-    B3 = 0.0
+    RealT = eltype(x)
+    x_0 = -4
+    rho = x[1] <= x_0 ? RealT(3.5) : 1 + RealT(0.2) * sin(5 * x[1])
+    v1 = x[1] <= x_0 ? RealT(5.8846) : zero(RealT)
+    v2 = x[1] <= x_0 ? RealT(1.1198) : zero(RealT)
+    v3 = 0
+    p = x[1] <= x_0 ? RealT(42.0267) : one(RealT)
+    B1 = 1
+    B2 = x[1] <= x_0 ? RealT(3.6359) : one(RealT)
+    B3 = 0
 
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end
@@ -45,15 +46,16 @@ function initial_condition_shu_osher_shock_tube_flipped(x, t,
                                                         equations::IdealGlmMhdEquations1D)
     # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
     # initial shock location is taken to be at x = 4
-    x_0 = 4.0
-    rho = x[1] <= x_0 ? 1.0 + 0.2 * sin(5.0 * x[1]) : 3.5
-    v1 = x[1] <= x_0 ? 0.0 : -5.8846
-    v2 = x[1] <= x_0 ? 0.0 : -1.1198
-    v3 = 0.0
-    p = x[1] <= x_0 ? 1.0 : 42.0267
-    B1 = 1.0
-    B2 = x[1] <= x_0 ? 1.0 : 3.6359
-    B3 = 0.0
+    RealT = eltype(x)
+    x_0 = 4
+    rho = x[1] <= x_0 ? 1 + RealT(0.2) * sin(5 * x[1]) : RealT(3.5)
+    v1 = x[1] <= x_0 ? zero(RealT) : RealT(-5.8846)
+    v2 = x[1] <= x_0 ? zero(RealT) : RealT(-1.1198)
+    v3 = 0
+    p = x[1] <= x_0 ? one(RealT) : RealT(42.0267)
+    B1 = 1
+    B2 = x[1] <= x_0 ? one(RealT) : RealT(3.6359)
+    B3 = 0
 
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end