Upwind FDSBP cleanup, part 2 (#1285)

* implementation of positive/negative part preferred by Michael Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * more experimental warnings * let Julia handle promotion of numbers * improve docstring of FluxUpwind * FDSBP constructor * use D_upw as canonical name in examples Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
trixi-framework · Dec 2, 2022 · 16385b3 · 16385b3
1 parent 545783c
commit 16385b3
Show file tree

Hide file tree

Showing 22 changed files with 130 additions and 90 deletions.
diff --git a/examples/tree_1d_fdsbp/elixir_burgers_basic.jl b/examples/tree_1d_fdsbp/elixir_burgers_basic.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -18,9 +17,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=32)
 flux_splitting = splitting_lax_friedrichs
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = 0.0
 coordinates_max = 1.0

diff --git a/examples/tree_1d_fdsbp/elixir_burgers_linear_stability.jl b/examples/tree_1d_fdsbp/elixir_burgers_linear_stability.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -21,9 +20,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_lax_friedrichs
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = -1.0
 coordinates_max =  1.0

diff --git a/examples/tree_1d_fdsbp/elixir_euler_convergence.jl b/examples/tree_1d_fdsbp/elixir_euler_convergence.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -18,9 +17,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                           xmin=-1.0, xmax=1.0,
                           N=32)
 flux_splitting = splitting_steger_warming
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = 0.0
 coordinates_max = 2.0

diff --git a/examples/tree_1d_fdsbp/elixir_euler_density_wave.jl b/examples/tree_1d_fdsbp/elixir_euler_density_wave.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -17,9 +16,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_coirier_vanleer
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = -1.0
 coordinates_max =  1.0

diff --git a/examples/tree_2d_fdsbp/elixir_advection_extended.jl b/examples/tree_2d_fdsbp/elixir_advection_extended.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -16,9 +15,9 @@ initial_condition = initial_condition_convergence_test
 D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
                             derivative_order=1, accuracy_order=4,
                             xmin=0.0, xmax=1.0, N=100)
-solver = DG(D_SBP, nothing #= mortar =#,
-            SurfaceIntegralStrongForm(flux_lax_friedrichs),
-            VolumeIntegralStrongForm())
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral=VolumeIntegralStrongForm())
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = ( 1.0,  1.0)

diff --git a/examples/tree_2d_fdsbp/elixir_euler_convergence.jl b/examples/tree_2d_fdsbp/elixir_euler_convergence.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -12,15 +11,15 @@ equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
 
-D = upwind_operators(SummationByPartsOperators.Mattsson2017,
-                     derivative_order=1,
-                     accuracy_order=4,
-                     xmin=-1.0, xmax=1.0,
-                     N=16)
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=16)
 flux_splitting = splitting_steger_warming
-solver = DG(D, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = ( 1.0,  1.0)

diff --git a/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -31,9 +30,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_vanleer_haenel
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = ( 1.0,  1.0)

diff --git a/examples/tree_2d_fdsbp/elixir_euler_vortex.jl b/examples/tree_2d_fdsbp/elixir_euler_vortex.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -61,9 +60,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_steger_warming
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = (-10.0, -10.0)
 coordinates_max = ( 10.0,  10.0)

diff --git a/examples/tree_3d_fdsbp/elixir_euler_convergence.jl b/examples/tree_3d_fdsbp/elixir_euler_convergence.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -18,9 +17,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_steger_warming
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = (0.0, 0.0, 0.0)
 coordinates_max = (2.0, 2.0, 2.0)

diff --git a/examples/tree_3d_fdsbp/elixir_euler_taylor_green_vortex.jl b/examples/tree_3d_fdsbp/elixir_euler_taylor_green_vortex.jl
@@ -1,6 +1,5 @@
-# TODO: FD
-# !!! warning "Experimental feature"
-#     This is an experimental feature and may change in any future releases.
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -36,9 +35,9 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin=-1.0, xmax=1.0,
                          N=16)
 flux_splitting = splitting_steger_warming
-solver = DG(D_upw, nothing #= mortar =#,
-            SurfaceIntegralUpwind(flux_splitting),
-            VolumeIntegralUpwind(flux_splitting))
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
 
 coordinates_min = (-1.0, -1.0, -1.0) .* pi
 coordinates_max = ( 1.0,  1.0,  1.0) .* pi

diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -189,6 +189,7 @@ export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,
+       FDSBP,
        VolumeIntegralWeakForm, VolumeIntegralStrongForm,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume,

diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl
@@ -199,8 +199,7 @@ 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)
+  return max(x, zero(x))
 end
 
 """
@@ -210,8 +209,7 @@ 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)
+  return min(x, zero(x))
 end
 
 

diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
@@ -374,6 +374,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - Joseph L. Steger and R. F. Warming (1979)
@@ -460,6 +463,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - Bram van Leer (1982)
@@ -551,6 +557,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - William Coirier and Bram van Leer (1991)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
@@ -308,11 +308,11 @@ Should be used together with [`UnstructuredMesh2D`](@ref).
   # [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
-    p_star = p_local * (1.0 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * equations.gamma * equations.inv_gamma_minus_one)
+    p_star = p_local * (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * equations.gamma * equations.inv_gamma_minus_one)
   else # v_normal > 0.0
-    A = 2.0 / ((equations.gamma + 1) * rho_local)
+    A = 2 / ((equations.gamma + 1) * rho_local)
     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.0 * A * (p_local + B)))
+    p_star = p_local + 0.5 * v_normal / A * (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
@@ -678,6 +678,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - Joseph L. Steger and R. F. Warming (1979)
@@ -804,6 +807,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - Bram van Leer (1982)
@@ -899,6 +905,9 @@ 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}`.
+
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
 """
 @inline function splitting_lax_friedrichs(u, orientation::Integer,
                                           equations::CompressibleEulerEquations2D)

diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
@@ -308,11 +308,11 @@ Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of
   # [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
-    p_star = p_local * (1.0 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * equations.gamma * equations.inv_gamma_minus_one)
+    p_star = p_local * (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * equations.gamma * equations.inv_gamma_minus_one)
   else # v_normal > 0.0
-    A = 2.0 / ((equations.gamma + 1) * rho_local)
+    A = 2 / ((equations.gamma + 1) * rho_local)
     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.0 * A * (p_local + B)))
+    p_star = p_local + 0.5 * v_normal / A * (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
@@ -694,6 +694,9 @@ 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}`.
 
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+
 ## References
 
 - Joseph L. Steger and R. F. Warming (1979)

diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl
@@ -145,6 +145,9 @@ 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}`.
+
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
 """
 @inline function splitting_lax_friedrichs(u, orientation::Integer,
                                           equations::InviscidBurgersEquation1D)

diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl
@@ -324,9 +324,12 @@ end
 A numerical flux `f(u_left, u_right) = f⁺(u_left) + f⁻(u_right)` based on
 flux vector splitting.
 
-The [`SurfaceIntegralUpwind`](@ref) with a given `splitting` is analytically
-equivalent to the [`SurfaceIntegralStrongForm`](@ref) with `FluxUpwind(splitting)`
-as numerical flux.
+The [`SurfaceIntegralUpwind`](@ref) with a given `splitting` is equivalent to
+the [`SurfaceIntegralStrongForm`](@ref) with `FluxUpwind(splitting)`
+as numerical flux (up to floating point differences).
+
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
 """
 struct FluxUpwind{Splitting}
   splitting::Splitting