diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl
index 7007bea887a..998deb04de2 100644
--- a/src/equations/shallow_water_1d.jl
+++ b/src/equations/shallow_water_1d.jl
@@ -71,16 +71,16 @@ A smooth initial condition used for convergence tests in combination with
 [`source_terms_convergence_test`](@ref)
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 """
-
 function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations1D)
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
-    c = 7.0
-    omega_x = 2.0 * pi * sqrt(2.0)
-    omega_t = 2.0 * pi
+    RealT = eltype(x)
+    c = 7
+    omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2))
+    omega_t = 2 * convert(RealT, pi)
 
     H = c + cos(omega_x * x[1]) * cos(omega_t * t)
-    v = 0.5
-    b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1])
+    v = 0.5f0
+    b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x[1])
     return prim2cons(SVector(H, v, b), equations)
 end
 
@@ -92,19 +92,19 @@ Source terms used for convergence tests in combination with
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 
 This manufactured solution source term is specifically designed for the bottom topography function
-`b(x) = 2.0 + 0.5 * sin(sqrt(2.0)*pi*x[1])`
+`b(x) = 2.0 + 0.5 * sinpi(sqrt(2.0) * x[1])`
 as defined in [`initial_condition_convergence_test`](@ref).
 """
-
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::ShallowWaterEquations1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocity is taken to be constant
-    c = 7.0
-    omega_x = 2.0 * pi * sqrt(2.0)
-    omega_t = 2.0 * pi
-    omega_b = sqrt(2.0) * pi
-    v = 0.5
+    RealT = eltype(u)
+    c = 7
+    omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2))
+    omega_t = 2 * convert(RealT, pi)
+    omega_b = sqrt(convert(RealT, 2)) * convert(RealT, pi)
+    v = 0.5f0
 
     sinX, cosX = sincos(omega_x * x[1])
     sinT, cosT = sincos(omega_t * t)
@@ -116,12 +116,12 @@ as defined in [`initial_condition_convergence_test`](@ref).
     H_t = -omega_t * cosX * sinT
 
     # bottom topography and its spatial derivative
-    b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1])
-    b_x = 0.5 * omega_b * cos(omega_b * x[1])
+    b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x[1])
+    b_x = 0.5f0 * omega_b * cos(omega_b * x[1])
 
     du1 = H_t + v * (H_x - b_x)
     du2 = v * du1 + equations.gravity * (H - b) * H_x
-    return SVector(du1, du2, 0.0)
+    return SVector(du1, du2, 0)
 end
 
 """
@@ -131,13 +131,14 @@ A weak blast wave discontinuity useful for testing, e.g., total energy conservat
 Note for the shallow water equations to the total energy acts as a mathematical entropy function.
 """
 function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations1D)
-    inicenter = 0.7
+    RealT = eltype(x)
+    inicenter = convert(RealT, 0.7)
     x_norm = x[1] - inicenter
     r = abs(x_norm)
 
     # Calculate primitive variables
-    H = r > 0.5 ? 3.25 : 4.0
-    v = r > 0.5 ? 0.0 : 0.1882
+    H = r > 0.5f0 ? 3.25f0 : 4.0f0
+    v = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882)
     b = sin(x[1]) # arbitrary continuous function
 
     return prim2cons(SVector(H, v, b), equations)
@@ -185,12 +186,12 @@ end
     h, h_v, _ = u
     v = velocity(u, equations)
 
-    p = 0.5 * equations.gravity * h^2
+    p = 0.5f0 * equations.gravity * h^2
 
     f1 = h_v
     f2 = h_v * v + p
 
-    return SVector(f1, f2, zero(eltype(u)))
+    return SVector(f1, f2, 0)
 end
 
 """
@@ -212,10 +213,8 @@ Further details are available in the paper:
     h_ll = waterheight(u_ll, equations)
     b_rr = u_rr[3]
 
-    z = zero(eltype(u_ll))
-
     # Bottom gradient nonconservative term: (0, g h b_x, 0)
-    f = SVector(z, equations.gravity * h_ll * b_rr, z)
+    f = SVector(0, equations.gravity * h_ll * b_rr, 0)
 
     return f
 end
@@ -249,19 +248,17 @@ and for curvilinear 2D case in the paper:
     h_ll, _, b_ll = u_ll
     h_rr, _, b_rr = u_rr
 
-    h_average = 0.5 * (h_ll + h_rr)
+    h_average = 0.5f0 * (h_ll + h_rr)
     b_jump = b_rr - b_ll
 
     # Includes two parts:
     #  (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
     #       cross-averaging across a discontinuous bottom topography
     #  (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
-    z = zero(eltype(u_ll))
-
-    f = SVector(z,
+    f = SVector(0,
                 equations.gravity * h_ll * b_ll +
                 equations.gravity * h_average * b_jump,
-                z)
+                0)
 
     return f
 end
@@ -296,15 +293,14 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    z = zero(eltype(u_ll))
     # Includes two parts:
     #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
     #        cross-averaging across a discontinuous bottom topography
     #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
-    return SVector(z,
+    return SVector(0,
                    equations.gravity * h_ll * b_ll +
                    equations.gravity * (h_ll^2 - h_ll_star^2),
-                   z)
+                   0)
 end
 
 """
@@ -337,10 +333,8 @@ For further details see:
     # Calculate jump
     b_jump = b_rr - b_ll
 
-    z = zero(eltype(u_ll))
-
     # Bottom gradient nonconservative term: (0, g h b_x, 0)
-    f = SVector(z, equations.gravity * h_ll * b_jump, z)
+    f = SVector(0, equations.gravity * h_ll * b_jump, 0)
 
     return f
 end
@@ -367,15 +361,15 @@ Details are available in Eq. (4.1) in the paper:
     v_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    h_avg = 0.5 * (h_ll + h_rr)
-    v_avg = 0.5 * (v_ll + v_rr)
-    p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2)
+    h_avg = 0.5f0 * (h_ll + h_rr)
+    v_avg = 0.5f0 * (v_ll + v_rr)
+    p_avg = 0.25f0 * equations.gravity * (h_ll^2 + h_rr^2)
 
     # Calculate fluxes depending on orientation
     f1 = h_avg * v_avg
     f2 = f1 * v_avg + p_avg
 
-    return SVector(f1, f2, zero(eltype(u_ll)))
+    return SVector(f1, f2, 0)
 end
 
 """
@@ -403,14 +397,14 @@ Further details are available in Theorem 1 of the paper:
     v_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    v_avg = 0.5 * (v_ll + v_rr)
-    p_avg = 0.5 * equations.gravity * h_ll * h_rr
+    v_avg = 0.5f0 * (v_ll + v_rr)
+    p_avg = 0.5f0 * equations.gravity * h_ll * h_rr
 
     # Calculate fluxes depending on orientation
-    f1 = 0.5 * (h_v_ll + h_v_rr)
+    f1 = 0.5f0 * (h_v_ll + h_v_rr)
     f2 = f1 * v_avg + p_avg
 
-    return SVector(f1, f2, zero(eltype(u_ll)))
+    return SVector(f1, f2, 0)
 end
 
 """
@@ -438,8 +432,8 @@ Further details on this hydrostatic reconstruction and its motivation can be fou
     v1_rr = velocity(u_rr, equations)
 
     # Compute the reconstructed water heights
-    h_ll_star = max(zero(h_ll), h_ll + b_ll - max(b_ll, b_rr))
-    h_rr_star = max(zero(h_rr), h_rr + b_rr - max(b_ll, b_rr))
+    h_ll_star = max(0, h_ll + b_ll - max(b_ll, b_rr))
+    h_rr_star = max(0, h_rr + b_rr - max(b_ll, b_rr))
 
     # Create the conservative variables using the reconstruted water heights
     u_ll_star = SVector(h_ll_star, h_ll_star * v1_ll, b_ll)
@@ -471,8 +465,8 @@ end
                                                               equations::ShallowWaterEquations1D)
     λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
                                   equations)
-    diss = -0.5 * λ * (u_rr - u_ll)
-    return SVector(diss[1], diss[2], zero(eltype(u_ll)))
+    diss = -0.5f0 * λ * (u_rr - u_ll)
+    return SVector(diss[1], diss[2], 0)
 end
 
 # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography
@@ -494,7 +488,7 @@ end
         factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min
         diss = u_rr - u_ll
         return factor_ll * f_ll - factor_rr * f_rr +
-               factor_diss * SVector(diss[1], diss[2], zero(eltype(u_ll)))
+               factor_diss * SVector(diss[1], diss[2], 0)
     end
 end
 
@@ -581,7 +575,7 @@ end
 
     v = velocity(u, equations)
 
-    w1 = equations.gravity * (h + b) - 0.5 * v^2
+    w1 = equations.gravity * (h + b) - 0.5f0 * v^2
     w2 = v
 
     return SVector(w1, w2, b)
@@ -591,7 +585,7 @@ end
 @inline function entropy2cons(w, equations::ShallowWaterEquations1D)
     w1, w2, b = w
 
-    h = (w1 + 0.5 * w2^2) / equations.gravity - b
+    h = (w1 + 0.5f0 * w2^2) / equations.gravity - b
     h_v = h * w2
     return SVector(h, h_v, b)
 end
@@ -612,7 +606,7 @@ end
 
 @inline function pressure(u, equations::ShallowWaterEquations1D)
     h = waterheight(u, equations)
-    p = 0.5 * equations.gravity * h^2
+    p = 0.5f0 * equations.gravity * h^2
     return p
 end
 
@@ -638,7 +632,7 @@ Or equation (9.17) in [this lecture notes](https://metaphor.ethz.ch/x/2019/hs/40
     h_rr = waterheight(u_rr, equations)
     v_rr = velocity(u_rr, equations)
 
-    h_roe = 0.5 * (h_ll + h_rr)
+    h_roe = 0.5f0 * (h_ll + h_rr)
     c_roe = sqrt(equations.gravity * h_roe)
 
     h_ll_sqrt = sqrt(h_ll)
@@ -658,7 +652,7 @@ end
 @inline function energy_total(cons, equations::ShallowWaterEquations1D)
     h, h_v, b = cons
 
-    e = (h_v^2) / (2 * h) + 0.5 * equations.gravity * h^2 + equations.gravity * h * b
+    e = (h_v^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b
     return e
 end
 
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index 73fae89a0fa..db8f00fc15b 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -77,16 +77,18 @@ A smooth initial condition used for convergence tests in combination with
 """
 function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations2D)
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2
-    c = 7.0
-    omega_x = 2.0 * pi * sqrt(2.0)
-    omega_t = 2.0 * pi
+    RealT = eltype(x)
+    c = 7
+    omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2))
+    omega_t = 2 * convert(RealT, pi)
 
     x1, x2 = x
 
     H = c + cos(omega_x * x1) * sin(omega_x * x2) * cos(omega_t * t)
-    v1 = 0.5
-    v2 = 1.5
-    b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2)
+    v1 = 0.5f0
+    v2 = 1.5f0
+    b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x1) +
+        0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x2)
     return prim2cons(SVector(H, v1, v2, b), equations)
 end
 
@@ -98,19 +100,20 @@ Source terms used for convergence tests in combination with
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 
 This manufactured solution source term is specifically designed for the bottom topography function
-`b(x,y) = 2 + 0.5 * sin(sqrt(2)*pi*x) + 0.5 * sin(sqrt(2)*pi*y)`
+`b(x,y) = 2 + 0.5 * sinpi(sqrt(2) * x) + 0.5 * sinpi(sqrt(2) * y)`
 as defined in [`initial_condition_convergence_test`](@ref).
 """
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::ShallowWaterEquations2D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocities are taken to be constants
-    c = 7.0
-    omega_x = 2.0 * pi * sqrt(2.0)
-    omega_t = 2.0 * pi
-    omega_b = sqrt(2.0) * pi
-    v1 = 0.5
-    v2 = 1.5
+    RealT = eltype(u)
+    c = 7
+    omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2))
+    omega_t = 2 * convert(RealT, pi)
+    omega_b = sqrt(convert(RealT, 2)) * convert(RealT, pi)
+    v1 = 0.5f0
+    v2 = 1.5f0
 
     x1, x2 = x
 
@@ -126,15 +129,16 @@ as defined in [`initial_condition_convergence_test`](@ref).
     H_t = -omega_t * cosX * sinY * sinT
 
     # bottom topography and its gradient
-    b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2)
-    tmp1 = 0.5 * omega_b
+    b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x1) +
+        0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x2)
+    tmp1 = 0.5f0 * omega_b
     b_x = tmp1 * cos(omega_b * x1)
     b_y = tmp1 * cos(omega_b * x2)
 
     du1 = H_t + v1 * (H_x - b_x) + v2 * (H_y - b_y)
     du2 = v1 * du1 + equations.gravity * (H - b) * H_x
     du3 = v2 * du1 + equations.gravity * (H - b) * H_y
-    return SVector(du1, du2, du3, 0.0)
+    return SVector(du1, du2, du3, 0)
 end
 
 """
@@ -145,7 +149,8 @@ Note for the shallow water equations to the total energy acts as a mathematical
 """
 function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations2D)
     # Set up polar coordinates
-    inicenter = SVector(0.7, 0.7)
+    RealT = eltype(x)
+    inicenter = SVector(convert(RealT, 0.7), convert(RealT, 0.7))
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
@@ -153,10 +158,10 @@ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquation
     sin_phi, cos_phi = sincos(phi)
 
     # Calculate primitive variables
-    H = r > 0.5 ? 3.25 : 4.0
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
-    b = 0.0 # by default assume there is no bottom topography
+    H = r > 0.5f0 ? 3.25f0 : 4.0f0
+    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi
+    b = 0 # by default assume there is no bottom topography
 
     return prim2cons(SVector(H, v1, v2, b), equations)
 end
@@ -185,8 +190,8 @@ For details see Section 9.2.5 of the book:
 
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
-                         u_inner[2] - 2.0 * u_normal * normal[1],
-                         u_inner[3] - 2.0 * u_normal * normal[2],
+                         u_inner[2] - 2 * u_normal * normal[1],
+                         u_inner[3] - 2 * u_normal * normal[2],
                          u_inner[4])
 
     # calculate the boundary flux
@@ -228,7 +233,7 @@ end
     h, h_v1, h_v2, _ = u
     v1, v2 = velocity(u, equations)
 
-    p = 0.5 * equations.gravity * h^2
+    p = 0.5f0 * equations.gravity * h^2
     if orientation == 1
         f1 = h_v1
         f2 = h_v1 * v1 + p
@@ -238,7 +243,7 @@ end
         f2 = h_v2 * v1
         f3 = h_v2 * v2 + p
     end
-    return SVector(f1, f2, f3, zero(eltype(u)))
+    return SVector(f1, f2, f3, 0)
 end
 
 # Calculate 1D flux for a single point in the normal direction
@@ -250,12 +255,12 @@ end
 
     v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
     h_v_normal = h * v_normal
-    p = 0.5 * equations.gravity * h^2
+    p = 0.5f0 * equations.gravity * h^2
 
     f1 = h_v_normal
     f2 = h_v_normal * v1 + p * normal_direction[1]
     f3 = h_v_normal * v2 + p * normal_direction[2]
-    return SVector(f1, f2, f3, zero(eltype(u)))
+    return SVector(f1, f2, f3, 0)
 end
 
 """
@@ -286,12 +291,11 @@ Further details are available in the paper:
     h_ll = waterheight(u_ll, equations)
     b_rr = u_rr[4]
 
-    z = zero(eltype(u_ll))
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     if orientation == 1
-        f = SVector(z, equations.gravity * h_ll * b_rr, z, z)
+        f = SVector(0, equations.gravity * h_ll * b_rr, 0, 0)
     else # orientation == 2
-        f = SVector(z, z, equations.gravity * h_ll * b_rr, z)
+        f = SVector(0, 0, equations.gravity * h_ll * b_rr, 0)
     end
     return f
 end
@@ -305,10 +309,10 @@ end
     b_rr = u_rr[4]
     # Note this routine only uses the `normal_direction_average` and the average of the
     # bottom topography to get a quadratic split form DG gradient on curved elements
-    return SVector(zero(eltype(u_ll)),
+    return SVector(0,
                    normal_direction_average[1] * equations.gravity * h_ll * b_rr,
                    normal_direction_average[2] * equations.gravity * h_ll * b_rr,
-                   zero(eltype(u_ll)))
+                   0)
 end
 
 """
@@ -349,24 +353,23 @@ and for curvilinear 2D case in the paper:
     h_ll, _, _, b_ll = u_ll
     h_rr, _, _, b_rr = u_rr
 
-    h_average = 0.5 * (h_ll + h_rr)
+    h_average = 0.5f0 * (h_ll + h_rr)
     b_jump = b_rr - b_ll
 
     # Includes two parts:
     #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
     #        cross-averaging across a discontinuous bottom topography
     #   (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
-    z = zero(eltype(u_ll))
     if orientation == 1
-        f = SVector(z,
+        f = SVector(0,
                     equations.gravity * h_ll * b_ll +
                     equations.gravity * h_average * b_jump,
-                    z, z)
+                    0, 0)
     else # orientation == 2
-        f = SVector(z, z,
+        f = SVector(0, 0,
                     equations.gravity * h_ll * b_ll +
                     equations.gravity * h_average * b_jump,
-                    z)
+                    0)
     end
 
     return f
@@ -389,14 +392,14 @@ end
 
     #   (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump`
     #        to handle discontinuous bathymetry
-    h_average = 0.5 * (h_ll + h_rr)
+    h_average = 0.5f0 * (h_ll + h_rr)
     b_jump = b_rr - b_ll
 
     f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump
     f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump
 
     # First and last equations do not have a nonconservative flux
-    f1 = f4 = zero(eltype(u_ll))
+    f1 = f4 = 0
 
     return SVector(f1, f2, f3, f4)
 end
@@ -426,8 +429,8 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
     v1_rr, v2_rr = velocity(u_rr, equations)
 
     # Compute the reconstructed water heights
-    h_ll_star = max(zero(h_ll), h_ll + b_ll - max(b_ll, b_rr))
-    h_rr_star = max(zero(h_rr), h_rr + b_rr - max(b_ll, b_rr))
+    h_ll_star = max(0, h_ll + b_ll - max(b_ll, b_rr))
+    h_rr_star = max(0, h_rr + b_rr - max(b_ll, b_rr))
 
     # Create the conservative variables using the reconstruted water heights
     u_ll_star = SVector(h_ll_star, h_ll_star * v1_ll, h_ll_star * v2_ll, b_ll)
@@ -469,21 +472,20 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    z = zero(eltype(u_ll))
     # Includes two parts:
     #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
     #        cross-averaging across a discontinuous bottom topography
     #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
     if orientation == 1
-        f = SVector(z,
+        f = SVector(0,
                     equations.gravity * h_ll * b_ll +
                     equations.gravity * (h_ll^2 - h_ll_star^2),
-                    z, z)
+                    0, 0)
     else # orientation == 2
-        f = SVector(z, z,
+        f = SVector(0, 0,
                     equations.gravity * h_ll * b_ll +
                     equations.gravity * (h_ll^2 - h_ll_star^2),
-                    z)
+                    0)
     end
 
     return f
@@ -516,7 +518,7 @@ end
     f3 += normal_direction_ll[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
 
     # First and last equations do not have a nonconservative flux
-    f1 = f4 = zero(eltype(u_ll))
+    f1 = f4 = 0
 
     return SVector(f1, f2, f3, f4)
 end
@@ -560,12 +562,11 @@ For further details see:
     # Calculate jump
     b_jump = b_rr - b_ll
 
-    z = zero(eltype(u_ll))
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     if orientation == 1
-        f = SVector(z, equations.gravity * h_ll * b_jump, z, z)
+        f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0)
     else # orientation == 2
-        f = SVector(z, z, equations.gravity * h_ll * b_jump, z)
+        f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0)
     end
     return f
 end
@@ -583,10 +584,10 @@ end
     b_jump = b_rr - b_ll
     # Note this routine only uses the `normal_direction_average` and the average of the
     # bottom topography to get a quadratic split form DG gradient on curved elements
-    return SVector(zero(eltype(u_ll)),
+    return SVector(0,
                    normal_direction_average[1] * equations.gravity * h_ll * b_jump,
                    normal_direction_average[2] * equations.gravity * h_ll * b_jump,
-                   zero(eltype(u_ll)))
+                   0)
 end
 
 """
@@ -611,10 +612,10 @@ Details are available in Eq. (4.1) in the paper:
     v1_rr, v2_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    h_avg = 0.5 * (h_ll + h_rr)
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2)
+    h_avg = 0.5f0 * (h_ll + h_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.25f0 * equations.gravity * (h_ll^2 + h_rr^2)
 
     # Calculate fluxes depending on orientation
     if orientation == 1
@@ -627,7 +628,7 @@ Details are available in Eq. (4.1) in the paper:
         f3 = f1 * v2_avg + p_avg
     end
 
-    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+    return SVector(f1, f2, f3, 0)
 end
 
 @inline function flux_fjordholm_etal(u_ll, u_rr, normal_direction::AbstractVector,
@@ -642,19 +643,19 @@ end
     v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
 
     # Average each factor of products in flux
-    h_avg = 0.5 * (h_ll + h_rr)
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    h2_avg = 0.5 * (h_ll^2 + h_rr^2)
-    p_avg = 0.5 * equations.gravity * h2_avg
-    v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
+    h_avg = 0.5f0 * (h_ll + h_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    h2_avg = 0.5f0 * (h_ll^2 + h_rr^2)
+    p_avg = 0.5f0 * equations.gravity * h2_avg
+    v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
 
     # Calculate fluxes depending on normal_direction
     f1 = h_avg * v_dot_n_avg
     f2 = f1 * v1_avg + p_avg * normal_direction[1]
     f3 = f1 * v2_avg + p_avg * normal_direction[2]
 
-    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+    return SVector(f1, f2, f3, 0)
 end
 
 """
@@ -682,22 +683,22 @@ Further details are available in Theorem 1 of the paper:
     v1_rr, v2_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.5 * equations.gravity * h_ll * h_rr
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * equations.gravity * h_ll * h_rr
 
     # Calculate fluxes depending on orientation
     if orientation == 1
-        f1 = 0.5 * (h_v1_ll + h_v1_rr)
+        f1 = 0.5f0 * (h_v1_ll + h_v1_rr)
         f2 = f1 * v1_avg + p_avg
         f3 = f1 * v2_avg
     else
-        f1 = 0.5 * (h_v2_ll + h_v2_rr)
+        f1 = 0.5f0 * (h_v2_ll + h_v2_rr)
         f2 = f1 * v1_avg
         f3 = f1 * v2_avg + p_avg
     end
 
-    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+    return SVector(f1, f2, f3, 0)
 end
 
 @inline function flux_wintermeyer_etal(u_ll, u_rr, normal_direction::AbstractVector,
@@ -711,18 +712,18 @@ end
     v1_rr, v2_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    h_v1_avg = 0.5 * (h_v1_ll + h_v1_rr)
-    h_v2_avg = 0.5 * (h_v2_ll + h_v2_rr)
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.5 * equations.gravity * h_ll * h_rr
+    h_v1_avg = 0.5f0 * (h_v1_ll + h_v1_rr)
+    h_v2_avg = 0.5f0 * (h_v2_ll + h_v2_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * equations.gravity * h_ll * h_rr
 
     # Calculate fluxes depending on normal_direction
     f1 = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2]
     f2 = f1 * v1_avg + p_avg * normal_direction[1]
     f3 = f1 * v2_avg + p_avg * normal_direction[2]
 
-    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+    return SVector(f1, f2, f3, 0)
 end
 
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
@@ -771,8 +772,8 @@ end
                                                               equations::ShallowWaterEquations2D)
     λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
                                   equations)
-    diss = -0.5 * λ * (u_rr - u_ll)
-    return SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll)))
+    diss = -0.5f0 * λ * (u_rr - u_ll)
+    return SVector(diss[1], diss[2], diss[3], 0)
 end
 
 # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography
@@ -794,7 +795,7 @@ end
         factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min
         diss = u_rr - u_ll
         return factor_ll * f_ll - factor_rr * f_rr +
-               factor_diss * SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll)))
+               factor_diss * SVector(diss[1], diss[2], diss[3], 0)
     end
 end
 
@@ -958,7 +959,7 @@ end
     v1, v2 = velocity(u, equations)
     v_square = v1^2 + v2^2
 
-    w1 = equations.gravity * (h + b) - 0.5 * v_square
+    w1 = equations.gravity * (h + b) - 0.5f0 * v_square
     w2 = v1
     w3 = v2
     return SVector(w1, w2, w3, b)
@@ -968,7 +969,7 @@ end
 @inline function entropy2cons(w, equations::ShallowWaterEquations2D)
     w1, w2, w3, b = w
 
-    h = (w1 + 0.5 * (w2^2 + w3^2)) / equations.gravity - b
+    h = (w1 + 0.5f0 * (w2^2 + w3^2)) / equations.gravity - b
     h_v1 = h * w2
     h_v2 = h * w3
     return SVector(h, h_v1, h_v2, b)
@@ -990,7 +991,7 @@ end
 
 @inline function pressure(u, equations::ShallowWaterEquations2D)
     h = waterheight(u, equations)
-    p = 0.5 * equations.gravity * h^2
+    p = 0.5f0 * equations.gravity * h^2
     return p
 end
 
@@ -1017,7 +1018,7 @@ slides 8 and 9.
     h_rr = waterheight(u_rr, equations)
     v1_rr, v2_rr = velocity(u_rr, equations)
 
-    h_roe = 0.5 * (h_ll + h_rr)
+    h_roe = 0.5f0 * (h_ll + h_rr)
     c_roe = sqrt(equations.gravity * h_roe)
 
     h_ll_sqrt = sqrt(h_ll)
@@ -1041,7 +1042,7 @@ end
 
     norm_ = norm(normal_direction)
 
-    h_roe = 0.5 * (h_ll + h_rr)
+    h_roe = 0.5f0 * (h_ll + h_rr)
     c_roe = sqrt(equations.gravity * h_roe) * norm_
 
     h_ll_sqrt = sqrt(h_ll)
@@ -1064,7 +1065,7 @@ end
 @inline function energy_total(cons, equations::ShallowWaterEquations2D)
     h, h_v1, h_v2, b = cons
 
-    e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5 * equations.gravity * h^2 +
+    e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 +
         equations.gravity * h * b
     return e
 end
diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 08620021254..eb90fb53f32 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -72,11 +72,12 @@ function initial_condition_convergence_test(x, t,
                                             equations::ShallowWaterEquationsQuasi1D)
     # generates a manufactured solution. 
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
-    Omega = sqrt(2) * pi
-    H = 2.0 + 0.5 * sin(Omega * x[1] - t)
-    v = 0.25
-    b = 0.2 - 0.05 * sin(Omega * x[1])
-    a = 1 + 0.1 * cos(Omega * x[1])
+    RealT = eltype(x)
+    Omega = sqrt(convert(RealT, 2)) * convert(RealT, pi)
+    H = 2 + 0.5f0 * sin(Omega * x[1] - t)
+    v = 0.25f0
+    b = convert(RealT, 0.2) - convert(RealT, 0.05) * sin(Omega * x[1])
+    a = 1 + convert(RealT, 0.1) * cos(Omega * x[1])
     return prim2cons(SVector(H, v, b, a), equations)
 end
 
@@ -88,30 +89,31 @@ Source terms used for convergence tests in combination with
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 
 This manufactured solution source term is specifically designed for the bottom topography function
-`b(x) = 0.2 - 0.05 * sin(sqrt(2) * pi *x[1])` and channel width 'a(x)= 1 + 0.1 * cos(sqrt(2) * pi * x[1])'
+`b(x) = 0.2 - 0.05 * sinpi(sqrt(2) * x[1])` and channel width 'a(x)= 1 + 0.1 * cospi(sqrt(2) * x[1])'
 as defined in [`initial_condition_convergence_test`](@ref).
 """
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::ShallowWaterEquationsQuasi1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocity is taken to be constant
-    Omega = sqrt(2) * pi
-    H = 2.0 + 0.5 * sin(Omega * x[1] - t)
-    H_x = 0.5 * cos(Omega * x[1] - t) * Omega
-    H_t = -0.5 * cos(Omega * x[1] - t)
+    RealT = eltype(u)
+    Omega = sqrt(convert(RealT, 2)) * convert(RealT, pi)
+    H = 2 + 0.5f0 * sin(Omega * x[1] - t)
+    H_x = 0.5f0 * cos(Omega * x[1] - t) * Omega
+    H_t = -0.5f0 * cos(Omega * x[1] - t)
 
-    v = 0.25
+    v = 0.25f0
 
-    b = 0.2 - 0.05 * sin(Omega * x[1])
-    b_x = -0.05 * cos(Omega * x[1]) * Omega
+    b = convert(RealT, 0.2) - convert(RealT, 0.05) * sin(Omega * x[1])
+    b_x = -convert(RealT, 0.05) * cos(Omega * x[1]) * Omega
 
-    a = 1 + 0.1 * cos(Omega * x[1])
-    a_x = -0.1 * sin(Omega * x[1]) * Omega
+    a = 1 + convert(RealT, 0.1) * cos(Omega * x[1])
+    a_x = -convert(RealT, 0.1) * sin(Omega * x[1]) * Omega
 
     du1 = a * H_t + v * (a_x * (H - b) + a * (H_x - b_x))
     du2 = v * du1 + a * (equations.gravity * (H - b) * H_x)
 
-    return SVector(du1, du2, 0.0, 0.0)
+    return SVector(du1, du2, 0, 0)
 end
 
 # Calculate 1D conservative flux for a single point
@@ -123,7 +125,7 @@ end
     f1 = a_h_v
     f2 = a_h_v * v
 
-    return SVector(f1, f2, zero(eltype(u)), zero(eltype(u)))
+    return SVector(f1, f2, 0, 0)
 end
 
 """
@@ -153,9 +155,7 @@ Further details are available in the paper:
     h_ll = waterheight(u_ll, equations)
     h_rr = waterheight(u_rr, equations)
 
-    z = zero(eltype(u_ll))
-
-    return SVector(z, equations.gravity * a_ll * h_ll * (h_rr + b_rr), z, z)
+    return SVector(0, equations.gravity * a_ll * h_ll * (h_rr + b_rr), 0, 0)
 end
 
 # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
@@ -200,10 +200,10 @@ Further details are available in the paper:
     v_ll = velocity(u_ll, equations)
     v_rr = velocity(u_rr, equations)
 
-    f1 = 0.5 * (a_h_v_ll + a_h_v_rr)
-    f2 = f1 * 0.5 * (v_ll + v_rr)
+    f1 = 0.5f0 * (a_h_v_ll + a_h_v_rr)
+    f2 = f1 * 0.5f0 * (v_ll + v_rr)
 
-    return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll)))
+    return SVector(f1, f2, 0, 0)
 end
 
 # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
@@ -240,8 +240,8 @@ end
                                                               equations::ShallowWaterEquationsQuasi1D)
     λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
                                   equations)
-    diss = -0.5 * λ * (u_rr - u_ll)
-    return SVector(diss[1], diss[2], zero(eltype(u_ll)), zero(eltype(u_ll)))
+    diss = -0.5f0 * λ * (u_rr - u_ll)
+    return SVector(diss[1], diss[2], 0, 0)
 end
 
 @inline function max_abs_speeds(u, equations::ShallowWaterEquationsQuasi1D)
@@ -278,7 +278,7 @@ end
     h = waterheight(u, equations)
     v = velocity(u, equations)
     #entropy variables are the same as ones in standard shallow water equations
-    w1 = equations.gravity * (h + b) - 0.5 * v^2
+    w1 = equations.gravity * (h + b) - 0.5f0 * v^2
     w2 = v
 
     return SVector(w1, w2, b, a)
@@ -306,7 +306,7 @@ end
 # Calculate total energy for a conservative state `cons`
 @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = cons
-    e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) +
+    e = (a_h_v^2) / (2 * a * a_h) + 0.5f0 * equations.gravity * (a_h^2 / a) +
         equations.gravity * a_h * b
     return e
 end
diff --git a/test/test_type.jl b/test/test_type.jl
index 24c005d1408..a3e50b81cfb 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -1852,6 +1852,320 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test typeof(@inferred energy_total(u, equations)) == RealT
         end
     end
+
+    @timed_testset "Shallow Water 1D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81))
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT))
+            orientation = 1
+            directions = [1, 2]
+            normal_direction = SVector(one(RealT))
+
+            surface_flux_function = flux_lax_friedrichs
+            dissipation = DissipationLocalLaxFriedrichs()
+            numflux = FluxHLL()
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
+                  RealT
+
+            for direction in directions
+                if RealT == Float32
+                    # check `surface_flux_function` (test broken)
+                    @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                               orientation,
+                                                                               direction,
+                                                                               x, t,
+                                                                               surface_flux_function,
+                                                                               equations)) ==
+                                 RealT
+                else
+                    @test eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                        orientation,
+                                                                        direction,
+                                                                        x, t,
+                                                                        surface_flux_function,
+                                                                        equations)) == RealT
+                end
+
+                @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, direction,
+                                                                equations)) ==
+                      RealT
+            end
+
+            @test eltype(@inferred flux(u, orientation, equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
+                                                                         orientation,
+                                                                         equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr,
+                                                                       orientation,
+                                                                       equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr,
+                                                                     orientation,
+                                                                     equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, orientation,
+                                                                    equations)) == RealT
+            @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation,
+                                                       equations)) == RealT
+            @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation,
+                                                         equations)) == RealT
+
+            @test eltype(eltype(@inferred hydrostatic_reconstruction_audusse_etal(u_ll,
+                                                                                  u_rr,
+                                                                                  equations))) ==
+                  RealT
+
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT
+            @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred numflux(u_ll, u_rr, orientation, equations)) == RealT
+            # no matching method 
+            # @test eltype(@inferred numflux(u_ll, u_rr, normal_direction, equations)) == RealT
+            @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation,
+                                                          equations)) == RealT
+            @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+
+            @test typeof(@inferred velocity(u, equations)) == RealT
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred prim2cons(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+            @test eltype(@inferred entropy2cons(u, equations)) == RealT
+            @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT
+            @test typeof(@inferred pressure(u, equations)) == RealT
+            @test typeof(@inferred waterheight_pressure(u, equations)) == RealT
+
+            @test typeof(@inferred entropy(cons, equations)) == RealT
+            @test typeof(@inferred energy_total(cons, equations)) == RealT
+            @test typeof(@inferred energy_kinetic(u, equations)) == RealT
+            @test typeof(@inferred energy_internal(cons, equations)) == RealT
+            @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT
+        end
+    end
+
+    @timed_testset "Shallow Water 2D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred ShallowWaterEquations2D(gravity_constant = RealT(9.81))
+
+            x = SVector(zero(RealT), zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT),
+                                                       one(RealT))
+            orientations = [1, 2]
+            directions = [1, 2, 3, 4]
+            normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT),
+                                                                                        zero(RealT))
+
+            surface_flux_function = flux_lax_friedrichs
+            dissipation = DissipationLocalLaxFriedrichs()
+            numflux = FluxHLL()
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
+                  RealT
+
+            if RealT == Float32
+                # check `surface_flux_function` (test broken)
+                @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                           normal_direction,
+                                                                           x, t,
+                                                                           surface_flux_function,
+                                                                           equations)) ==
+                             RealT
+            else
+                @test eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                    normal_direction,
+                                                                    x, t,
+                                                                    surface_flux_function,
+                                                                    equations)) == RealT
+            end
+            @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
+                                                                         normal_direction_ll,
+                                                                         normal_direction_average,
+                                                                         equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr,
+                                                                       normal_direction_ll,
+                                                                       normal_direction_average,
+                                                                       equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr,
+                                                                     normal_direction_ll,
+                                                                     normal_direction_average,
+                                                                     equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                                    normal_direction_ll,
+                                                                    normal_direction_average,
+                                                                    equations)) == RealT
+            @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, normal_direction,
+                                                         equations)) == RealT
+
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) ==
+                  RealT
+            @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred numflux(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction,
+                                                          equations)) == RealT
+            @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, normal_direction,
+                                                            equations)) == RealT
+
+            for orientation in orientations
+                for direction in directions
+                    if RealT == Float32
+                        # check `surface_flux_function` (test broken)
+                        @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                                   orientation,
+                                                                                   direction,
+                                                                                   x, t,
+                                                                                   surface_flux_function,
+                                                                                   equations)) ==
+                                     RealT
+                    else
+                        @test eltype(@inferred boundary_condition_slip_wall(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_nonconservative_wintermeyer_etal(u_ll, u_rr,
+                                                                             orientation,
+                                                                             equations)) ==
+                      RealT
+                @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr,
+                                                                           orientation,
+                                                                           equations)) ==
+                      RealT
+                @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr,
+                                                                         orientation,
+                                                                         equations)) ==
+                      RealT
+                @test eltype(eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                                               orientation,
+                                                                               equations))) ==
+                      RealT
+                @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation,
+                                                           equations)) == RealT
+                @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation,
+                                                             equations)) == RealT
+
+                @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation,
+                                                           equations)) ==
+                      RealT
+                @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) ==
+                      RealT
+                @test eltype(@inferred numflux(u_ll, u_rr, orientation, equations)) == RealT
+                @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation,
+                                                           equations)) == RealT
+                @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation,
+                                                           equations)) == RealT
+                @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation,
+                                                              equations)) == RealT
+                @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, orientation,
+                                                                equations)) == RealT
+            end
+
+            @test eltype(eltype(@inferred hydrostatic_reconstruction_audusse_etal(u_ll,
+                                                                                  u_rr,
+                                                                                  equations))) ==
+                  RealT
+
+            @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+            @test eltype(@inferred velocity(u, equations)) == RealT
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred prim2cons(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+            @test eltype(@inferred entropy2cons(u, equations)) == RealT
+            @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT
+            @test typeof(@inferred pressure(u, equations)) == RealT
+            @test typeof(@inferred waterheight_pressure(u, equations)) == RealT
+
+            @test typeof(@inferred entropy(cons, equations)) == RealT
+            @test typeof(@inferred energy_total(cons, equations)) == RealT
+            @test typeof(@inferred energy_kinetic(u, equations)) == RealT
+            @test typeof(@inferred energy_internal(cons, equations)) == RealT
+            @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT
+        end
+    end
+
+    @timed_testset "Shallow Water Quasi 1D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred ShallowWaterEquationsQuasi1D(gravity_constant = RealT(9.81))
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT))
+            orientation = 1
+            normal_direction = normal_ll = normal_rr = SVector(one(RealT))
+
+            dissipation = DissipationLocalLaxFriedrichs()
+
+            @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                  RealT
+            @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
+                  RealT
+
+            @test eltype(@inferred flux(u, orientation, equations)) == RealT
+            @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                                  orientation,
+                                                                  equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                                  normal_direction,
+                                                                  equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll,
+                                                                  normal_rr,
+                                                                  equations)) == RealT
+            @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation,
+                                                  equations)) == RealT
+            @test eltype(@inferred flux_chan_etal(u_ll, u_rr, normal_direction,
+                                                  equations)) == RealT
+
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) ==
+                  RealT
+            @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT
+            @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+            @test typeof(@inferred velocity(u, equations)) == RealT
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred prim2cons(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+            @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT
+
+            @test typeof(@inferred entropy(cons, equations)) == RealT
+            @test typeof(@inferred energy_total(cons, equations)) == RealT
+            @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT
+        end
+    end
 end
 
 end # module