Upwind FDSBP performance improvements (#1282)

* improve some docstrings (#1274)

* fix convergence_test for elixirs without trailing new line (#1280)

* Fix `initial_condition_isentropic_vortex` (#1279)

* strengthen the isentropic vortex initial conditions in TreeMesh elixirs

* update test values for new isentropic vortex initial conditions

* update vortex initial condition in special elixir and docs

* fix typos in the IC

* update tree_mpi test values

* remove comment lines because them seem to break literate

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* improve performance of 3D volume integral with steger_warming_splitting by ca. 1/4 for Taylor-Green (serial)

* lose 2 percent of threaded volume terms performance but simplify code and reduce duplications

* remove flux_upwind stub

* 2D and 1D as well; improvement ca. 20% in 3D, 10% in 2D, not much in 1D

* positive_part, negative_part for 3D Steger-Warming; improves serial PID of TGV by ca. 10%

* positive_part, negative_part for 2D, 1D

* simplify CSE for LLVM (ca. 2 % for TGV)

* comments on f_minus_plus etc.

* improve comment

* comments on positive/negative part

Co-authored-by: Andrew Winters <[email protected]>
Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

src/auxiliary/math.jl

@@ -192,5 +192,27 @@ julia> min(2, 5, 1)
 
 
 
+"""
+    positive_part(x)
+
+Return `x` if `x` is positive, else zero. In other words, return
+`(x + abs(x)) / 2` for real numbers `x`.
+"""
+@inline function positive_part(x)
+  o = zero(x)
+  return ifelse(x > o, x, o)
+end
+
+"""
+    negative_part(x)
+
+Return `x` if `x` is negative, else zero. In other words, return
+`(x - abs(x)) / 2` for real numbers `x`.
+"""
+@inline function negative_part(x)
+  o = zero(x)
+  return ifelse(x < o, x, o)
+end
+
 
 end # @muladd

src/equations/compressible_euler_1d.jl

@@ -361,16 +361,18 @@ end
 
 
 """
+    splitting_steger_warming(u, orientation::Integer,
+                             equations::CompressibleEulerEquations1D)
     splitting_steger_warming(u, which::Union{Val{:minus}, Val{:plus}}
                              orientation::Integer,
                              equations::CompressibleEulerEquations1D)
 
 Splitting of the compressible Euler flux of Steger and Warming.
 
-Returns the flux "minus" (associated with waves going into the
-negative axis direction) or "plus" (associated with waves going into the
-positive axis direction), determined by the argument `which` set to
-`Val{:minus}()` or `Val{:plus}`.
+Returns a tuple of the fluxes "minus" (associated with waves going into the
+negative axis direction) and "plus" (associated with waves going into the
+positive axis direction). If only one of the fluxes is required, use the
+function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`.
 
 ## References
 
@@ -379,6 +381,13 @@ positive axis direction), determined by the argument `which` set to
   With Application to Finite Difference Methods
   [NASA Technical Memorandum](https://ntrs.nasa.gov/api/citations/19790020779/downloads/19790020779.pdf)
 """
+@inline function splitting_steger_warming(u, orientation::Integer,
+                                          equations::CompressibleEulerEquations1D)
+  fm = splitting_steger_warming(u, Val{:minus}(), orientation, equations)
+  fp = splitting_steger_warming(u, Val{:plus}(),  orientation, equations)
+  return fm, fp
+end
+
 @inline function splitting_steger_warming(u, ::Val{:plus}, orientation::Integer,
                                           equations::CompressibleEulerEquations1D)
   rho, rho_v1, rho_e = u
@@ -390,16 +399,17 @@ positive axis direction), determined by the argument `which` set to
   lambda2 = v1 + a
   lambda3 = v1 - a
 
-  lambda1_p = 0.5 * (lambda1 + abs(lambda1))
-  lambda2_p = 0.5 * (lambda2 + abs(lambda2))
-  lambda3_p = 0.5 * (lambda3 + abs(lambda3))
+  lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :)
+  lambda2_p = positive_part(lambda2)
+  lambda3_p = positive_part(lambda3)
 
   alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p
 
-  f1p = 0.5 * rho / equations.gamma * alpha_p
-  f2p = 0.5 * rho / equations.gamma * (alpha_p * v1 + a * (lambda2_p - lambda3_p))
-  f3p = 0.5 * rho / equations.gamma * (alpha_p * 0.5 * v1^2 + a * v1 * (lambda2_p - lambda3_p)
-                                       + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one)
+  rho_2gamma = 0.5 * rho / equations.gamma
+  f1p = rho_2gamma * alpha_p
+  f2p = rho_2gamma * (alpha_p * v1 + a * (lambda2_p - lambda3_p))
+  f3p = rho_2gamma * (alpha_p * 0.5 * v1^2 + a * v1 * (lambda2_p - lambda3_p)
+                      + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one)
 
   return SVector(f1p, f2p, f3p)
 end
@@ -415,22 +425,25 @@ end
     lambda2 = v1 + a
     lambda3 = v1 - a
 
-    lambda1_m = 0.5 * (lambda1 - abs(lambda1))
-    lambda2_m = 0.5 * (lambda2 - abs(lambda2))
-    lambda3_m = 0.5 * (lambda3 - abs(lambda3))
+    lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :)
+    lambda2_m = negative_part(lambda2)
+    lambda3_m = negative_part(lambda3)
 
     alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m
 
-    f1m = 0.5 * rho / equations.gamma * alpha_m
-    f2m = 0.5 * rho / equations.gamma * (alpha_m * v1 + a * (lambda2_m - lambda3_m))
-    f3m = 0.5 * rho / equations.gamma * (alpha_m * 0.5 * v1^2 + a * v1 * (lambda2_m - lambda3_m)
-                                         + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one)
+    rho_2gamma = 0.5 * rho / equations.gamma
+    f1m = rho_2gamma * alpha_m
+    f2m = rho_2gamma * (alpha_m * v1 + a * (lambda2_m - lambda3_m))
+    f3m = rho_2gamma * (alpha_m * 0.5 * v1^2 + a * v1 * (lambda2_m - lambda3_m)
+                        + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one)
 
   return SVector(f1m, f2m, f3m)
 end
 
 
 """
+    splitting_vanleer_haenel(u, orientation::Integer,
+                             equations::CompressibleEulerEquations1D)
     splitting_vanleer_haenel(u, which::Union{Val{:minus}, Val{:plus}}
                              orientation::Integer,
                              equations::CompressibleEulerEquations1D)
@@ -442,10 +455,10 @@ convective terms. As such there are many pressure splittings suggested across
 the literature. We implement the 'p4' variant suggested by Liou and Steffen as
 it proved the most robust in practice.
 
-Returns the flux "minus" (associated with waves going into the
-negative axis direction) or "plus" (associated with waves going into the
-positive axis direction), determined by the argument `which` set to
-`Val{:minus}()` or `Val{:plus}`.
+Returns a tuple of the fluxes "minus" (associated with waves going into the
+negative axis direction) and "plus" (associated with waves going into the
+positive axis direction). If only one of the fluxes is required, use the
+function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`.
 
 ## References
 
@@ -459,6 +472,13 @@ positive axis direction), determined by the argument `which` set to
   High-Order Polynomial Expansions (HOPE) for Flux-Vector Splitting
   [NASA Technical Memorandum](https://ntrs.nasa.gov/citations/19910016425)
 """
+@inline function splitting_vanleer_haenel(u, orientation::Integer,
+                                          equations::CompressibleEulerEquations1D)
+  fm = splitting_vanleer_haenel(u, Val{:minus}(), orientation, equations)
+  fp = splitting_vanleer_haenel(u, Val{:plus}(),  orientation, equations)
+  return fm, fp
+end
+
 @inline function splitting_vanleer_haenel(u, ::Val{:plus}, orientation::Integer,
                                           equations::CompressibleEulerEquations1D)
   rho, rho_v1, rho_e = u
@@ -515,6 +535,8 @@ end
 # of this splitting on "el diablo" still crashes early. Can we learn anything
 # from the design of this splitting?
 """
+    splitting_coirier_vanleer(u, orientation::Integer,
+                              equations::CompressibleEulerEquations1D)
     splitting_coirier_vanleer(u, which::Union{Val{:minus}, Val{:plus}}
                               orientation::Integer,
                               equations::CompressibleEulerEquations1D)
@@ -524,10 +546,10 @@ The splitting has correction terms in the pressure splitting as well as
 the mass and energy flux components. The motivation for these corrections
 are to handle flows at the low Mach number limit.
 
-Returns the flux "minus" (associated with waves going into the
-negative axis direction) or "plus" (associated with waves going into the
-positive axis direction), determined by the argument `which` set to
-`Val{:minus}()` or `Val{:plus}`.
+Returns a tuple of the fluxes "minus" (associated with waves going into the
+negative axis direction) and "plus" (associated with waves going into the
+positive axis direction). If only one of the fluxes is required, use the
+function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`.
 
 ## References
 
@@ -536,6 +558,13 @@ positive axis direction), determined by the argument `which` set to
   II - Progress in flux-vector splitting
   [DOI: 10.2514/6.1991-1566](https://doi.org/10.2514/6.1991-1566)
 """
+@inline function splitting_coirier_vanleer(u, orientation::Integer,
+                                          equations::CompressibleEulerEquations1D)
+  fm = splitting_coirier_vanleer(u, Val{:minus}(), orientation, equations)
+  fp = splitting_coirier_vanleer(u, Val{:plus}(),  orientation, equations)
+  return fm, fp
+end
+
 @inline function splitting_coirier_vanleer(u, ::Val{:plus}, orientation::Integer,
                                             equations::CompressibleEulerEquations1D)
   rho, rho_v1, rho_e = u