Upwind SBP on curved meshes

examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl

@@ -33,7 +33,7 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
                          xmin = -1.0, xmax = 1.0,
                          N = 9)
 
-flux_splitting = splitting_steger_warming
+flux_splitting = splitting_drikakis_tsangaris
 solver = FDSBP(D_upw,
                surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
                volume_integral = VolumeIntegralUpwind(flux_splitting))

src/Trixi.jl

@@ -191,7 +191,8 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        FluxUpwind
 
 export splitting_steger_warming, splitting_vanleer_haenel,
-       splitting_coirier_vanleer, splitting_lax_friedrichs
+       splitting_coirier_vanleer, splitting_lax_friedrichs,
+       splitting_drikakis_tsangaris
 
 export initial_condition_constant,
        initial_condition_gauss,

src/equations/compressible_euler_2d.jl

@@ -683,15 +683,15 @@ end
 end
 
 """
-    splitting_steger_warming(u, orientation_or_normal_direction,
+    splitting_steger_warming(u, orientation::Integer,
                              equations::CompressibleEulerEquations2D)
     splitting_steger_warming(u, which::Union{Val{:minus}, Val{:plus}}
-                             orientation_or_normal_direction,
+                             orientation::Integer,
                              equations::CompressibleEulerEquations2D)
 
 Splitting of the compressible Euler flux of Steger and Warming. For
-curvilinear coordinates to compute in the `normal_direction` we use
-a modified version of this splitting proposed by Drikakis and Tsangris.
+curvilinear coordinates use the improved Steger-Warming-type splitting
+[`splitting_drikakis_tsangaris`](@ref).
 
 Returns a tuple of the fluxes "minus" (associated with waves going into the
 negative axis direction) and "plus" (associated with waves going into the
@@ -707,12 +707,8 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()
   Flux Vector Splitting of the Inviscid Gasdynamic Equations
   With Application to Finite Difference Methods
   [NASA Technical Memorandum](https://ntrs.nasa.gov/api/citations/19790020779/downloads/19790020779.pdf)
-- D. Drikakis and S. Tsangaris (1993)
-  On the solution of the compressible Navier-Stokes equations using
-  improved flux vector splitting methods
-  [DOI: 10.1016/0307-904X(93)90054-K](https://doi.org/10.1016/0307-904X(93)90054-K)
 """
-@inline function splitting_steger_warming(u, orientation_or_normal_direction,
+@inline function splitting_steger_warming(u, orientation::Integer,
                                           equations::CompressibleEulerEquations2D)
     fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction,
                                   equations)
@@ -817,9 +813,122 @@ end
     return SVector(f1m, f2m, f3m, f4m)
 end
 
-@inline function splitting_steger_warming(u, ::Val{:plus},
-                                          normal_direction::AbstractVector,
-                                          equations::CompressibleEulerEquations2D)
+"""
+    splitting_drikakis_tsangaris(u, orientation_or_normal_direction,
+                                 equations::CompressibleEulerEquations2D)
+    splitting_drikakis_tsangaris(u, which::Union{Val{:minus}, Val{:plus}}
+                                 orientation_or_normal_direction,
+                                 equations::CompressibleEulerEquations2D)
+
+Improved variant of the Steger-Warming flux vector splitting for
+generalized coordinates. This splitting also reformulates the energy
+flux as in Hänel et al. to obtain conservation of the total temperature
+for inviscid flows.
+
+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.
+
+## References
+
+- D. Drikakis and S. Tsangaris (1993)
+  On the solution of the compressible Navier-Stokes equations using
+  improved flux vector splitting methods
+  [DOI: 10.1016/0307-904X(93)90054-K](https://doi.org/10.1016/0307-904X(93)90054-K)
+- D. Hänel, R. Schwane and G. Seider (1987)
+  On the accuracy of upwind schemes for the solution of the Navier-Stokes equations
+  [DOI: 10.2514/6.1987-1105](https://doi.org/10.2514/6.1987-1105)
+"""
+@inline function splitting_drikakis_tsangaris(u, orientation_or_normal_direction,
+                                              equations::CompressibleEulerEquations2D)
+    fm = splitting_drikakis_tsangaris(u, Val{:minus}(), orientation_or_normal_direction,
+                                      equations)
+    fp = splitting_drikakis_tsangaris(u, Val{:plus}(), orientation_or_normal_direction,
+                                      equations)
+    return fm, fp
+end
+
+@inline function splitting_drikakis_tsangaris(u, ::Val{:plus}, orientation::Integer,
+                                              equations::CompressibleEulerEquations2D)
+    rho, rho_v1, rho_v2, rho_e = u
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+    a = sqrt(equations.gamma * p / rho)
+    H = (rho_e + p) / rho
+
+    if orientation == 1
+        lambda1 = v1 + a
+        lambda2 = v1 - a
+
+        lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :)
+        lambda2_p = positive_part(lambda2)
+
+        rhoa_2gamma = 0.5 * rho * a / equations.gamma
+        f1p = 0.5 * rho * (lambda1_p + lambda2_p)
+        f2p = f1p * v1 + rhoa_2gamma * (lambda1_p - lambda2_p)
+        f3p = f1p * v2
+        f4p = f1p * H
+    else # orientation == 2
+        lambda1 = v2 + a
+        lambda2 = v2 - a
+
+        lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :)
+        lambda2_p = positive_part(lambda2)
+
+        rhoa_2gamma = 0.5 * rho * a / equations.gamma
+        f1p = 0.5 * rho * (lambda1_p + lambda2_p)
+        f2p = f1p * v1
+        f3p = f1p * v2 + rhoa_2gamma * (lambda1_p - lambda2_p)
+        f4p = f1p * H
+    end
+    return SVector(f1p, f2p, f3p, f4p)
+end
+
+@inline function splitting_drikakis_tsangaris(u, ::Val{:minus}, orientation::Integer,
+                                              equations::CompressibleEulerEquations2D)
+    rho, rho_v1, rho_v2, rho_e = u
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+    a = sqrt(equations.gamma * p / rho)
+    H = (rho_e + p) / rho
+
+    if orientation == 1
+        lambda1 = v1 + a
+        lambda2 = v1 - a
+
+        lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :)
+        lambda2_m = negative_part(lambda2)
+
+        rhoa_2gamma = 0.5 * rho * a / equations.gamma
+        f1m = 0.5 * rho * (lambda1_m + lambda2_m)
+        f2m = f1m * v1 + rhoa_2gamma * (lambda1_m - lambda2_m)
+        f3m = f1m * v2
+        f4m = f1m * H
+    else # orientation == 2
+        lambda1 = v2 + a
+        lambda2 = v2 - a
+
+        lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :)
+        lambda2_m = negative_part(lambda2)
+
+        rhoa_2gamma = 0.5 * rho * a / equations.gamma
+        f1m = 0.5 * rho * (lambda1_m + lambda2_m)
+        f2m = f1m * v1
+        f3m = f1m * v2 + rhoa_2gamma * (lambda1_m - lambda2_m)
+        f4m = f1m * H
+    end
+    return SVector(f1m, f2m, f3m, f4m)
+end
+
+@inline function splitting_drikakis_tsangaris(u, ::Val{:plus},
+                                              normal_direction::AbstractVector,
+                                              equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     v1 = rho_v1 / rho
     v2 = rho_v2 / rho
@@ -844,9 +953,9 @@ end
     return SVector(f1p, f2p, f3p, f4p)
 end
 
-@inline function splitting_steger_warming(u, ::Val{:minus},
-                                          normal_direction::AbstractVector,
-                                          equations::CompressibleEulerEquations2D)
+@inline function splitting_drikakis_tsangaris(u, ::Val{:minus},
+                                              normal_direction::AbstractVector,
+                                              equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     v1 = rho_v1 / rho
     v2 = rho_v2 / rho
@@ -975,7 +1084,8 @@ contains a reformulation due to Hänel et al. where the energy flux uses the
 enthalpy. The pressure splitting is independent from the splitting of the
 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.
+it proved the most robust in practice. For details on the curvilinear variant
+of this flux vector splitting see Anderson et al.
 
 Returns a tuple of the fluxes "minus" (associated with waves going into the
 negative axis direction) and "plus" (associated with waves going into the
@@ -996,6 +1106,9 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()
 - Meng-Sing Liou and Chris J. Steffen, Jr. (1991)
   High-Order Polynomial Expansions (HOPE) for Flux-Vector Splitting
   [NASA Technical Memorandum](https://ntrs.nasa.gov/citations/19910016425)
+- W. Kyle Anderson, James L. Thomas, and Bram van Leer (1986)
+  Comparison of Finite Volume Flux Vector Splittings for the Euler Equations
+  [DOI: 10.2514/3.9465](https://doi.org/10.2514/3.9465)
 """
 @inline function splitting_vanleer_haenel(u, orientation_or_normal_direction,
                                           equations::CompressibleEulerEquations2D)

src/solvers/fdsbp_unstructured/containers_2d.jl

@@ -9,21 +9,14 @@
 #! format: noindent
 
 # initialize all the values in the container of a general FD block (either straight sided or curved)
-# OBS! Requires the SBP derivative matrix in order to compute metric terms. For a standard
-# SBP operator the matrix is passed directly, whereas for an upwind SBP operator we must
-# specifically pass the central differencing matrix.
+# OBS! Requires the SBP derivative matrix in order to compute metric terms.
 function init_element!(elements, element, basis::AbstractDerivativeOperator,
                        corners_or_surface_curves)
     calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis),
                            corners_or_surface_curves)
 
-    if basis isa DerivativeOperator
-        calc_metric_terms!(elements.jacobian_matrix, element, basis,
-                           elements.node_coordinates)
-    else # basis isa SummationByPartsOperators.UpwindOperators
-        calc_metric_terms!(elements.jacobian_matrix, element, basis.central,
-                           elements.node_coordinates)
-    end
+    calc_metric_terms!(elements.jacobian_matrix, element, basis,
+                       elements.node_coordinates)
 
     calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
 
@@ -36,6 +29,13 @@ function init_element!(elements, element, basis::AbstractDerivativeOperator,
     return elements
 end
 
+# Specialization to pass the central differencing matrix from an upwind SBP operator
+function calc_metric_terms!(jacobian_matrix, element,
+                            D_SBP::SummationByPartsOperators.UpwindOperators,
+                            node_coordinates)
+    calc_metric_terms!(jacobian_matrix, element, D_SBP.central, node_coordinates)
+end
+
 # construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
 # applied to the node coordinates.
 function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator,

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

@@ -187,7 +187,7 @@ end
 # surface contributions are added.
 function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
                                 equations, surface_integral::SurfaceIntegralStrongForm,
-                                dg::FDSBP, cache)
+                                dg::DG, cache)
     inv_weight_left = inv(left_boundary_weight(dg.basis))
     inv_weight_right = inv(right_boundary_weight(dg.basis))
     @unpack normal_directions, surface_flux_values = cache.elements

test/test_tree_2d_fdsbp.jl

@@ -102,6 +102,32 @@ end
         end
     end
 
+    @trixi_testset "elixir_euler_convergence.jl with Drikakis-Tsangaris splitting" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence.jl"),
+                            l2=[
+                                1.708838999643608e-6,
+                                1.7437997854485807e-6,
+                                1.7437997854741082e-6,
+                                5.457223460116349e-6,
+                            ],
+                            linf=[
+                                9.796504911285808e-6,
+                                9.614745899888533e-6,
+                                9.614745899444443e-6,
+                                4.02610718399643e-5,
+                            ],
+                            tspan=(0.0, 0.1), flux_splitting=splitting_drikakis_tsangaris)
+
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
+
     @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_euler_kelvin_helmholtz_instability.jl"),