Own sqrt and log returning NaN for "correct" multi-thread behaviour

src/auxiliary/math.jl

@@ -5,6 +5,10 @@
 @muladd begin
 #! format: noindent
 
+sqrt_(x::T) where {T <: AbstractFloat} = x < zero(x) ? oftype(x, NaN) : Base.sqrt(x)
+# < 0 suffices since log(0) = -Inf
+log_(x::T) where {T <: AbstractFloat} = x < zero(x) ? oftype(x, NaN) : Base.log(x)
+
 """
     ln_mean(x, y)
 
@@ -59,7 +63,7 @@ Given ε = 1.0e-4, we use the following algorithm.
     if f2 < epsilon_f2
         return (x + y) / @evalpoly(f2, 2, 2/3, 2/5, 2/7)
     else
-        return (y - x) / log(y / x)
+        return (y - x) / log_(y / x)
     end
 end
 
@@ -79,7 +83,7 @@ multiplication.
     if f2 < epsilon_f2
         return @evalpoly(f2, 2, 2/3, 2/5, 2/7) / (x + y)
     else
-        return log(y / x) / (y - x)
+        return log_(y / x) / (y - x)
     end
 end
 

src/callbacks_step/averaging_dg2d.jl

@@ -68,8 +68,8 @@ function calc_mean_values!(mean_values, averaging_callback_cache, u, u_prev,
             rho, v1, v2, p = u_node_prim
             rho_prev, v1_prev, v2_prev, p_prev = u_prev_node_prim
 
-            c = sqrt(equations.gamma * p / rho)
-            c_prev = sqrt(equations.gamma * p_prev / rho_prev)
+            c = sqrt_(equations.gamma * p / rho)
+            c_prev = sqrt_(equations.gamma * p_prev / rho_prev)
 
             # Calculate the contribution to the mean values using the trapezoidal rule
             vorticity_mean[i, j, element] += integration_constant *

src/equations/compressible_euler_1d.jl

@@ -230,7 +230,7 @@ are available in the paper:
     # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
     # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
     if v_normal <= 0.0
-        sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed
+        sound_speed = sqrt_(equations.gamma * p_local / rho_local) # local sound speed
         p_star = p_local *
                  (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 *
                                                                              equations.gamma *
@@ -240,7 +240,7 @@ are available in the paper:
         B = p_local * (equations.gamma - 1) / (equations.gamma + 1)
         p_star = p_local +
                  0.5 * v_normal / A *
-                 (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B)))
+                 (v_normal + sqrt_(v_normal^2 + 4 * A * (p_local + B)))
     end
 
     # For the slip wall we directly set the flux as the normal velocity is zero
@@ -449,7 +449,7 @@ end
     rho, rho_v1, rho_e = u
     v1 = rho_v1 / rho
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
 
     lambda1 = v1
     lambda2 = v1 + a
@@ -475,7 +475,7 @@ end
     rho, rho_v1, rho_e = u
     v1 = rho_v1 / rho
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
 
     lambda1 = v1
     lambda2 = v1 + a
@@ -544,7 +544,7 @@ end
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
 
     # sound speed and enthalpy
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
     H = (rho_e + p) / rho
 
     # signed Mach number
@@ -566,7 +566,7 @@ end
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
 
     # sound speed and enthalpy
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
     H = (rho_e + p) / rho
 
     # signed Mach number
@@ -632,7 +632,7 @@ end
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
 
     # sound speed and enthalpy
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
     H = (rho_e + p) / rho
 
     # signed Mach number
@@ -659,7 +659,7 @@ end
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1)
 
     # sound speed and enthalpy
-    a = sqrt(equations.gamma * p / rho)
+    a = sqrt_(equations.gamma * p / rho)
     H = (rho_e + p) / rho
 
     # signed Mach number
@@ -690,11 +690,11 @@ end
     v1_ll = rho_v1_ll / rho_ll
     v_mag_ll = abs(v1_ll)
     p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v_mag_ll^2)
-    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_ll = sqrt_(equations.gamma * p_ll / rho_ll)
     v1_rr = rho_v1_rr / rho_rr
     v_mag_rr = abs(v1_rr)
     p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v_mag_rr^2)
-    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    c_rr = sqrt_(equations.gamma * p_rr / rho_rr)
 
     λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
 end
@@ -705,8 +705,8 @@ end
     rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
 
-    λ_min = v1_ll - sqrt(equations.gamma * p_ll / rho_ll)
-    λ_max = v1_rr + sqrt(equations.gamma * p_rr / rho_rr)
+    λ_min = v1_ll - sqrt_(equations.gamma * p_ll / rho_ll)
+    λ_max = v1_rr + sqrt_(equations.gamma * p_rr / rho_rr)
 
     return λ_min, λ_max
 end
@@ -717,8 +717,8 @@ end
     rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
 
-    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
-    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    c_ll = sqrt_(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt_(equations.gamma * p_rr / rho_rr)
 
     λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
     λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
@@ -742,20 +742,20 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     v1_ll = rho_v1_ll / rho_ll
     e_ll = rho_e_ll / rho_ll
     p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v1_ll^2)
-    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_ll = sqrt_(equations.gamma * p_ll / rho_ll)
 
     v1_rr = rho_v1_rr / rho_rr
     e_rr = rho_e_rr / rho_rr
     p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v1_rr^2)
-    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    c_rr = sqrt_(equations.gamma * p_rr / rho_rr)
 
     # Obtain left and right fluxes
     f_ll = flux(u_ll, orientation, equations)
     f_rr = flux(u_rr, orientation, equations)
 
     # Compute Roe averages
-    sqrt_rho_ll = sqrt(rho_ll)
-    sqrt_rho_rr = sqrt(rho_rr)
+    sqrt_rho_ll = sqrt_(rho_ll)
+    sqrt_rho_rr = sqrt_(rho_rr)
     sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr
     vel_L = v1_ll
     vel_R = v1_rr
@@ -764,7 +764,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     H_ll = (rho_e_ll + p_ll) / rho_ll
     H_rr = (rho_e_rr + p_rr) / rho_rr
     H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
-    c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
+    c_roe = sqrt_((equations.gamma - 1) * (H_roe - ekin_roe))
 
     Ssl = min(vel_L - c_ll, vel_roe - c_roe)
     Ssr = max(vel_R + c_rr, vel_roe + c_roe)
@@ -833,22 +833,22 @@ Compactly summarized:
 
     # `u_ll[3]` is total energy `rho_e_ll` on the left
     H_ll = (u_ll[3] + p_ll) / rho_ll
-    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_ll = sqrt_(equations.gamma * p_ll / rho_ll)
 
     # `u_rr[3]` is total energy `rho_e_rr` on the right
     H_rr = (u_rr[3] + p_rr) / rho_rr
-    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    c_rr = sqrt_(equations.gamma * p_rr / rho_rr)
 
     # Compute Roe averages
-    sqrt_rho_ll = sqrt(rho_ll)
-    sqrt_rho_rr = sqrt(rho_rr)
+    sqrt_rho_ll = sqrt_(rho_ll)
+    sqrt_rho_rr = sqrt_(rho_rr)
     inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr)
 
     v_roe = (sqrt_rho_ll * v_ll + sqrt_rho_rr * v_rr) * inv_sum_sqrt_rho
     v_roe_mag = v_roe^2
 
     H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho
-    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag))
+    c_roe = sqrt_((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag))
 
     # Compute convenience constant for positivity preservation, see
     # https://doi.org/10.1016/0021-9991(91)90211-3
@@ -865,7 +865,7 @@ end
     rho, rho_v1, rho_e = u
     v1 = rho_v1 / rho
     p = (equations.gamma - 1) * (rho_e - 1 / 2 * rho * v1^2)
-    c = sqrt(equations.gamma * p / rho)
+    c = sqrt_(equations.gamma * p / rho)
 
     return (abs(v1) + c,)
 end
@@ -887,7 +887,7 @@ end
     v1 = rho_v1 / rho
     v_square = v1^2
     p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square)
-    s = log(p) - equations.gamma * log(rho)
+    s = log_(p) - equations.gamma * log_(rho)
     rho_p = rho / p
 
     w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square
@@ -951,7 +951,7 @@ end
     p = (equations.gamma - 1) * (cons[3] - 1 / 2 * (cons[2]^2) / cons[1])
 
     # Thermodynamic entropy
-    s = log(p) - equations.gamma * log(cons[1])
+    s = log_(p) - equations.gamma * log_(cons[1])
 
     return s
 end