Add numerical support of other real types (compressible_euler)

src/equations/compressible_euler_multicomponent_1d.jl

@@ -123,14 +123,15 @@ A smooth initial condition used for convergence tests in combination with
 """
 function initial_condition_convergence_test(x, t,
                                             equations::CompressibleEulerMulticomponentEquations1D)
+    RealT = eltype(x)
     c = 2
-    A = 0.1
+    A = convert(RealT, 0.1)
     L = 2
-    f = 1 / L
-    omega = 2 * pi * f
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
     ini = c + A * sin(omega * (x[1] - t))
 
-    v1 = 1.0
+    v1 = 1
 
     rho = ini
 
@@ -159,20 +160,21 @@ Source terms used for convergence tests in combination with
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::CompressibleEulerMulticomponentEquations1D)
     # Same settings as in `initial_condition`
+    RealT = eltype(u)
     c = 2
-    A = 0.1
+    A = convert(RealT, 0.1)
     L = 2
-    f = 1 / L
-    omega = 2 * pi * f
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
 
     gamma = totalgamma(u, equations)
 
     x1, = x
     si, co = sincos((t - x1) * omega)
-    tmp = (-((4 * si * A - 4c) + 1) * (gamma - 1) * co * A * omega) / 2
+    tmp = (-((4 * si * A - 4 * c) + 1) * (gamma - 1) * co * A * omega) / 2
 
     # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1
-    du_rho = SVector{ncomponents(equations), real(equations)}(0.0
+    du_rho = SVector{ncomponents(equations), real(equations)}(0
                                                               for i in eachcomponent(equations))
 
     du1 = tmp
@@ -194,24 +196,25 @@ A for multicomponent adapted weak blast wave adapted to multicomponent and taken
 function initial_condition_weak_blast_wave(x, t,
                                            equations::CompressibleEulerMulticomponentEquations1D)
     # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
-    cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)
+    cos_phi = x_norm > 0 ? 1 : -1
 
-    prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ?
+    prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ?
                                                                 2^(i - 1) * (1 - 2) /
                                                                 (1 -
                                                                  2^ncomponents(equations)) *
-                                                                1.0 :
+                                                                1 :
                                                                 2^(i - 1) * (1 - 2) /
                                                                 (1 -
                                                                  2^ncomponents(equations)) *
-                                                                1.1691
+                                                                convert(RealT, 1.1691)
                                                                 for i in eachcomponent(equations))
 
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    p = r > 0.5 ? 1.0 : 1.245
+    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
+    p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245)
 
     prim_other = SVector{2, real(equations)}(v1, p)
 
@@ -227,7 +230,7 @@ end
 
     v1 = rho_v1 / rho
     gamma = totalgamma(u, equations)
-    p = (gamma - 1) * (rho_e - 0.5 * rho * v1^2)
+    p = (gamma - 1) * (rho_e - 0.5f0 * rho * v1^2)
 
     f_rho = densities(u, v1, equations)
     f1 = rho_v1 * v1 + p
@@ -255,7 +258,7 @@ Entropy conserving two-point flux by
     rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2],
                                                                          u_rr[i + 2])
                                                                  for i in eachcomponent(equations))
-    rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] +
+    rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] +
                                                                  u_rr[i + 2])
                                                                 for i in eachcomponent(equations))
 
@@ -269,28 +272,28 @@ Entropy conserving two-point flux by
     # extract velocities
     v1_ll = rho_v1_ll / rho_ll
     v1_rr = rho_v1_rr / rho_rr
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v1_square = 0.5 * (v1_ll^2 + v1_rr^2)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2)
     v_sum = v1_avg
 
-    enth = zero(v_sum)
-    help1_ll = zero(v1_ll)
-    help1_rr = zero(v1_rr)
+    enth = 0
+    help1_ll = 0
+    help1_rr = 0
 
     for i in eachcomponent(equations)
         enth += rhok_avg[i] * gas_constants[i]
         help1_ll += u_ll[i + 2] * cv[i]
         help1_rr += u_rr[i + 2] * cv[i]
     end
 
-    T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll
-    T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2)) / help1_rr
-    T = 0.5 * (1.0 / T_ll + 1.0 / T_rr)
-    T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr)
+    T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2)) / help1_ll
+    T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2)) / help1_rr
+    T = 0.5f0 * (1 / T_ll + 1 / T_rr)
+    T_log = ln_mean(1 / T_ll, 1 / T_rr)
 
     # Calculate fluxes depending on orientation
-    help1 = zero(T_ll)
-    help2 = zero(T_rr)
+    help1 = 0
+    help2 = 0
 
     f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg
                                                              for i in eachcomponent(equations))
@@ -299,7 +302,7 @@ Entropy conserving two-point flux by
         help2 += f_rho[i]
     end
     f1 = (help2) * v1_avg + enth / T
-    f2 = (help1) / T_log - 0.5 * (v1_square) * (help2) + v1_avg * f1
+    f2 = (help1) / T_log - 0.5f0 * (v1_square) * (help2) + v1_avg * f1
 
     f_other = SVector{2, real(equations)}(f1, f2)
 
@@ -330,7 +333,7 @@ See also
     rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2],
                                                                          u_rr[i + 2])
                                                                  for i in eachcomponent(equations))
-    rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] +
+    rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] +
                                                                  u_rr[i + 2])
                                                                 for i in eachcomponent(equations))
 
@@ -339,25 +342,25 @@ See also
     rho_rr = density(u_rr, equations)
 
     # Calculating gamma
-    gamma = totalgamma(0.5 * (u_ll + u_rr), equations)
+    gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations)
     inv_gamma_minus_one = 1 / (gamma - 1)
 
     # extract velocities
     v1_ll = rho_v1_ll / rho_ll
     v1_rr = rho_v1_rr / rho_rr
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    velocity_square_avg = 0.5 * (v1_ll * v1_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    velocity_square_avg = 0.5f0 * (v1_ll * v1_rr)
 
     # density flux
     f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg
                                                              for i in eachcomponent(equations))
 
     # helpful variables
-    f_rho_sum = zero(v1_ll)
-    help1_ll = zero(v1_ll)
-    help1_rr = zero(v1_rr)
-    enth_ll = zero(v1_ll)
-    enth_rr = zero(v1_rr)
+    f_rho_sum = 0
+    help1_ll = 0
+    help1_rr = 0
+    enth_ll = 0
+    enth_rr = 0
     for i in eachcomponent(equations)
         enth_ll += u_ll[i + 2] * gas_constants[i]
         enth_rr += u_rr[i + 2] * gas_constants[i]
@@ -367,17 +370,17 @@ See also
     end
 
     # temperature and pressure
-    T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll
-    T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2)) / help1_rr
+    T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2)) / help1_ll
+    T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2)) / help1_rr
     p_ll = T_ll * enth_ll
     p_rr = T_rr * enth_rr
-    p_avg = 0.5 * (p_ll + p_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
     inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
 
     # momentum and energy flux
     f1 = f_rho_sum * v1_avg + p_avg
     f2 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) +
-         0.5 * (p_ll * v1_rr + p_rr * v1_ll)
+         0.5f0 * (p_ll * v1_rr + p_rr * v1_ll)
     f_other = SVector{2, real(equations)}(f1, f2)
 
     return vcat(f_other, f_rho)
@@ -398,8 +401,8 @@ end
     v_ll = rho_v1_ll / rho_ll
     v_rr = rho_v1_rr / rho_rr
 
-    p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * rho_ll * v_ll^2)
-    p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * rho_rr * v_rr^2)
+    p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_ll^2)
+    p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * rho_rr * v_rr^2)
     c_ll = sqrt(gamma_ll * p_ll / rho_ll)
     c_rr = sqrt(gamma_rr * p_rr / rho_rr)
 
@@ -414,7 +417,7 @@ end
     v1 = rho_v1 / rho
 
     gamma = totalgamma(u, equations)
-    p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2))
+    p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2))
     c = sqrt(gamma * p / rho)
 
     return (abs(v1) + c,)
@@ -431,7 +434,7 @@ end
     v1 = rho_v1 / rho
     gamma = totalgamma(u, equations)
 
-    p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2))
+    p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2))
     prim_other = SVector{2, real(equations)}(v1, p)
 
     return vcat(prim_other, prim_rho)
@@ -451,7 +454,7 @@ end
 
     rho_v1 = rho * v1
 
-    rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1)
+    rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1)
 
     cons_other = SVector{2, RealT}(rho_v1, rho_e)
 
@@ -466,8 +469,8 @@ end
 
     rho = density(u, equations)
 
-    help1 = zero(rho)
-    gas_constant = zero(rho)
+    help1 = 0
+    gas_constant = 0
     for i in eachcomponent(equations)
         help1 += u[i + 2] * cv[i]
         gas_constant += gas_constants[i] * (u[i + 2] / rho)
@@ -477,10 +480,10 @@ end
     v_square = v1^2
     gamma = totalgamma(u, equations)
 
-    p = (gamma - 1) * (rho_e - 0.5 * rho * v_square)
+    p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square)
     s = log(p) - gamma * log(rho) - log(gas_constant)
     rho_p = rho / p
-    T = (rho_e - 0.5 * rho * v_square) / (help1)
+    T = (rho_e - 0.5f0 * rho * v_square) / (help1)
 
     entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] *
                                                                    (1 - log(T)) +
@@ -507,14 +510,14 @@ end
                                                                     (-cv[i] *
                                                                      log(-w[2]) -
                                                                      cp[i] + w[i + 2] -
-                                                                     0.5 * w[1]^2 /
+                                                                     0.5f0 * w[1]^2 /
                                                                      w[2]))
                                                                 for i in eachcomponent(equations))
 
-    rho = zero(cons_rho[1])
-    help1 = zero(cons_rho[1])
-    help2 = zero(cons_rho[1])
-    p = zero(cons_rho[1])
+    rho = 0
+    help1 = 0
+    help2 = 0
+    p = 0
     for i in eachcomponent(equations)
         rho += cons_rho[i]
         help1 += cons_rho[i] * cv[i] * gammas[i]
@@ -523,7 +526,7 @@ end
     end
     u1 = rho * v1
     gamma = help1 / help2
-    u2 = p / (gamma - 1) + 0.5 * rho * v1^2
+    u2 = p / (gamma - 1) + 0.5f0 * rho * v1^2
     cons_other = SVector{2, real(equations)}(u1, u2)
     return vcat(cons_other, cons_rho)
 end
@@ -534,7 +537,7 @@ end
     rho = density(u, equations)
     T = temperature(u, equations)
 
-    total_entropy = zero(u[1])
+    total_entropy = 0
     for i in eachcomponent(equations)
         total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2]))
     end
@@ -548,15 +551,15 @@ end
     rho_v1, rho_e = u
 
     rho = density(u, equations)
-    help1 = zero(rho)
+    help1 = 0
 
     for i in eachcomponent(equations)
         help1 += u[i + 2] * cv[i]
     end
 
     v1 = rho_v1 / rho
     v_square = v1^2
-    T = (rho_e - 0.5 * rho * v_square) / help1
+    T = (rho_e - 0.5f0 * rho * v_square) / help1
 
     return T
 end
@@ -570,8 +573,8 @@ partial density fractions as well as the partial specific heats at constant volu
 @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D)
     @unpack cv, gammas = equations
 
-    help1 = zero(u[1])
-    help2 = zero(u[1])
+    help1 = 0
+    help2 = 0
 
     for i in eachcomponent(equations)
         help1 += u[i + 2] * cv[i] * gammas[i]
@@ -587,13 +590,13 @@ end
     rho = density(u, equations)
     gamma = totalgamma(u, equations)
 
-    p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho)
+    p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho)
 
     return p
 end
 
 @inline function density(u, equations::CompressibleEulerMulticomponentEquations1D)
-    rho = zero(u[1])
+    rho = 0
 
     for i in eachcomponent(equations)
         rho += u[i + 2]