Implement upwind flux for linearized Euler equations (trixi-framework…

…#1557)

* Enable input checks for LEE keyword constructor

* Extend LEE implementation to curved meshes

* Implement upwind flux for linearized Euler equations

* Add upwind flux examples and tests

* Fix comments in linearized Euler elixirs

* Clarify LEE Gaussian source elixir

* Rename `flux_upwind` to `flux_godunov`

* Add parentheses around multiline expressions

* Add consistency checks for LEE Godunov flux

* Explain odd mean values in more detail

* Use normalized normal vector to simplify flux

* Add docstring for LEE upwind flux

* Update examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl

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

---------

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

examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl

@@ -0,0 +1,89 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+# Based on the TreeMesh example `elixir_acoustics_gaussian_source.jl`.
+# The acoustic perturbation equations have been replaced with the linearized Euler
+# equations and instead of the Cartesian `TreeMesh` a rotated `P4estMesh` is used
+
+# Oscillating Gaussian-shaped source terms
+function source_terms_gauss(u, x, t, equations::LinearizedEulerEquations2D)
+  r = 0.1
+  A = 1.0
+  f = 2.0
+
+  # Velocity sources
+  s2 = 0.0
+  s3 = 0.0
+  # Density and pressure source
+  s1 = s4 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sin(2 * pi * f * t)
+
+  return SVector(s1, s2, s3, s4)
+end
+
+initial_condition_zero(x, t, equations::LinearizedEulerEquations2D) = SVector(0.0, 0.0, 0.0, 0.0)
+
+###############################################################################
+# semidiscretization of the linearized Euler equations
+
+# Create a domain that is a 30° rotated version of [-3, 3]^2
+c = cospi(2 * 30.0 / 360.0)
+s = sinpi(2 * 30.0 / 360.0)
+rot_mat = Trixi.SMatrix{2, 2}([c -s; s c])
+mapping(xi, eta) = rot_mat * SVector(3.0*xi, 3.0*eta)
+
+# Mean density and speed of sound are slightly off from 1.0 to allow proper verification of
+# curved LEE implementation using this elixir (some things in the LEE cancel if both are 1.0)
+equations = LinearizedEulerEquations2D(v_mean_global=Tuple(rot_mat * SVector(-0.5, 0.25)),
+                                       c_mean_global=1.02, rho_mean_global=1.01)
+
+initial_condition = initial_condition_zero
+
+# Create DG solver with polynomial degree = 3 and upwind flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_godunov)
+
+# Create a uniformly refined mesh with periodic boundaries
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension, polydeg=1,
+                 mapping=mapping,
+                 periodicity=true, initial_refinement_level=2)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_gauss)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 2.0
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()

examples/tree_2d_dgsem/elixir_linearizedeuler_gauss_wall.jl

@@ -0,0 +1,68 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linearized Euler equations
+
+equations = LinearizedEulerEquations2D(v_mean_global=(0.5, 0.0), c_mean_global=1.0,
+                                       rho_mean_global=1.0)
+
+# Create DG solver with polynomial degree = 5 and upwind flux as surface flux
+solver = DGSEM(polydeg=5, surface_flux=flux_godunov)
+
+coordinates_min = (-100.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (100.0, 200.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=100_000,
+                periodicity=false)
+
+function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations2D)
+  v1_prime = 0.0
+  v2_prime = 0.0
+  rho_prime = p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25)
+  return SVector(rho_prime, v1_prime, v2_prime, p_prime)
+end
+initial_condition = initial_condition_gauss_wall
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_condition_wall)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 30.0
+tspan = (0.0, 30.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.7)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks)
+
+# Print the timer summary
+summary_callback()

src/equations/linearized_euler_2d.jl

@@ -53,7 +53,7 @@ end
 
 function LinearizedEulerEquations2D(; v_mean_global::NTuple{2, <:Real},
                                     c_mean_global::Real, rho_mean_global::Real)
-    return LinearizedEulerEquations2D(SVector(v_mean_global), c_mean_global,
+    return LinearizedEulerEquations2D(v_mean_global, c_mean_global,
                                       rho_mean_global)
 end
 
@@ -126,6 +126,24 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
+# Calculate 1D flux for a single point
+@inline function flux(u, normal_direction::AbstractVector,
+                      equations::LinearizedEulerEquations2D)
+    @unpack v_mean_global, c_mean_global, rho_mean_global = equations
+    rho_prime, v1_prime, v2_prime, p_prime = u
+
+    v_mean_normal = v_mean_global[1] * normal_direction[1] +
+                    v_mean_global[2] * normal_direction[2]
+    v_prime_normal = v1_prime * normal_direction[1] + v2_prime * normal_direction[2]
+
+    f1 = v_mean_normal * rho_prime + rho_mean_global * v_prime_normal
+    f2 = v_mean_normal * v1_prime + normal_direction[1] * p_prime / rho_mean_global
+    f3 = v_mean_normal * v2_prime + normal_direction[2] * p_prime / rho_mean_global
+    f4 = v_mean_normal * p_prime + c_mean_global^2 * rho_mean_global * v_prime_normal
+
+    return SVector(f1, f2, f3, f4)
+end
+
 @inline have_constant_speed(::LinearizedEulerEquations2D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations2D)
@@ -143,6 +161,198 @@ end
     end
 end
 
+@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::LinearizedEulerEquations2D)
+    @unpack v_mean_global, c_mean_global = equations
+    v_mean_normal = normal_direction[1] * v_mean_global[1] +
+                    normal_direction[2] * v_mean_global[2]
+    return abs(v_mean_normal) + c_mean_global * norm(normal_direction)
+end
+
+@doc raw"""
+    flux_godunov(u_ll, u_rr, orientation_or_normal_direction,
+                 equations::LinearizedEulerEquations2D)
+
+An upwind flux for the linearized Euler equations based on diagonalization of the physical
+flux matrix. Given the physical flux ``Au``, ``A=T \Lambda T^{-1}`` with
+``\Lambda`` being a diagonal matrix that holds the eigenvalues of ``A``, decompose
+``\Lambda = \Lambda^+ + \Lambda^-`` where ``\Lambda^+`` and ``\Lambda^-`` are diagonal
+matrices holding the positive and negative eigenvalues of ``A``, respectively. Then for
+left and right states ``u_L, u_R``, the numerical flux calculated by this function is given
+by ``A^+ u_L + A^- u_R`` where ``A^{\pm} = T \Lambda^{\pm} T^{-1}``.
+
+The diagonalization of the flux matrix can be found in
+- R. F. Warming, Richard M. Beam and B. J. Hyett (1975)
+  Diagonalization and simultaneous symmetrization of the gas-dynamic matrices
+  [DOI: 10.1090/S0025-5718-1975-0388967-5](https://doi.org/10.1090/S0025-5718-1975-0388967-5)
+"""
+@inline function flux_godunov(u_ll, u_rr, orientation::Integer,
+                              equations::LinearizedEulerEquations2D)
+    @unpack v_mean_global, rho_mean_global, c_mean_global = equations
+    v1_mean = v_mean_global[1]
+    v2_mean = v_mean_global[2]
+
+    rho_prime_ll, v1_prime_ll, v2_prime_ll, p_prime_ll = u_ll
+    rho_prime_rr, v1_prime_rr, v2_prime_rr, p_prime_rr = u_rr
+
+    if orientation == 1
+        # Eigenvalues of the flux matrix
+        lambda1 = v1_mean
+        lambda2 = v1_mean - c_mean_global
+        lambda3 = v1_mean + c_mean_global
+
+        lambda1_p = positive_part(lambda1)
+        lambda2_p = positive_part(lambda2)
+        lambda3_p = positive_part(lambda3)
+        lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
+        lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+
+        lambda1_m = negative_part(lambda1)
+        lambda2_m = negative_part(lambda2)
+        lambda3_m = negative_part(lambda3)
+        lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
+        lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+
+        f1p = (lambda1_p * rho_prime_ll +
+               lambda3m2_half_p / c_mean_global * rho_mean_global * v1_prime_ll +
+               (lambda2p3_half_p - lambda1_p) / c_mean_global^2 * p_prime_ll)
+        f2p = (lambda2p3_half_p * v1_prime_ll +
+               lambda3m2_half_p / c_mean_global * p_prime_ll / rho_mean_global)
+        f3p = lambda1_p * v2_prime_ll
+        f4p = (lambda3m2_half_p * c_mean_global * rho_mean_global * v1_prime_ll +
+               lambda2p3_half_p * p_prime_ll)
+
+        f1m = (lambda1_m * rho_prime_rr +
+               lambda3m2_half_m / c_mean_global * rho_mean_global * v1_prime_rr +
+               (lambda2p3_half_m - lambda1_m) / c_mean_global^2 * p_prime_rr)
+        f2m = (lambda2p3_half_m * v1_prime_rr +
+               lambda3m2_half_m / c_mean_global * p_prime_rr / rho_mean_global)
+        f3m = lambda1_m * v2_prime_rr
+        f4m = (lambda3m2_half_m * c_mean_global * rho_mean_global * v1_prime_rr +
+               lambda2p3_half_m * p_prime_rr)
+
+        f1 = f1p + f1m
+        f2 = f2p + f2m
+        f3 = f3p + f3m
+        f4 = f4p + f4m
+    else # orientation == 2
+        # Eigenvalues of the flux matrix
+        lambda1 = v2_mean
+        lambda2 = v2_mean - c_mean_global
+        lambda3 = v2_mean + c_mean_global
+
+        lambda1_p = positive_part(lambda1)
+        lambda2_p = positive_part(lambda2)
+        lambda3_p = positive_part(lambda3)
+        lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
+        lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+
+        lambda1_m = negative_part(lambda1)
+        lambda2_m = negative_part(lambda2)
+        lambda3_m = negative_part(lambda3)
+        lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
+        lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+
+        f1p = (lambda1_p * rho_prime_ll +
+               lambda3m2_half_p / c_mean_global * rho_mean_global * v2_prime_ll +
+               (lambda2p3_half_p - lambda1_p) / c_mean_global^2 * p_prime_ll)
+        f2p = lambda1_p * v1_prime_ll
+        f3p = (lambda2p3_half_p * v2_prime_ll +
+               lambda3m2_half_p / c_mean_global * p_prime_ll / rho_mean_global)
+        f4p = (lambda3m2_half_p * c_mean_global * rho_mean_global * v2_prime_ll +
+               lambda2p3_half_p * p_prime_ll)
+
+        f1m = (lambda1_m * rho_prime_rr +
+               lambda3m2_half_m / c_mean_global * rho_mean_global * v2_prime_rr +
+               (lambda2p3_half_m - lambda1_m) / c_mean_global^2 * p_prime_rr)
+        f2m = lambda1_m * v1_prime_rr
+        f3m = (lambda2p3_half_m * v2_prime_rr +
+               lambda3m2_half_m / c_mean_global * p_prime_rr / rho_mean_global)
+        f4m = (lambda3m2_half_m * c_mean_global * rho_mean_global * v2_prime_rr +
+               lambda2p3_half_m * p_prime_rr)
+
+        f1 = f1p + f1m
+        f2 = f2p + f2m
+        f3 = f3p + f3m
+        f4 = f4p + f4m
+    end
+
+    return SVector(f1, f2, f3, f4)
+end
+
+@inline function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
+                              equations::LinearizedEulerEquations2D)
+    @unpack v_mean_global, rho_mean_global, c_mean_global = equations
+    rho_prime_ll, v1_prime_ll, v2_prime_ll, p_prime_ll = u_ll
+    rho_prime_rr, v1_prime_rr, v2_prime_rr, p_prime_rr = u_rr
+
+    # Do not use `normalize` since we use `norm_` later to scale the eigenvalues
+    norm_ = norm(normal_direction)
+    normal_vector = normal_direction / norm_
+
+    # Use normalized vector here, scaling is applied via eigenvalues of the flux matrix
+    v_mean_normal = v_mean_global[1] * normal_vector[1] +
+                    v_mean_global[2] * normal_vector[2]
+    v_prime_normal_ll = v1_prime_ll * normal_vector[1] + v2_prime_ll * normal_vector[2]
+    v_prime_normal_rr = v1_prime_rr * normal_vector[1] + v2_prime_rr * normal_vector[2]
+
+    # Eigenvalues of the flux matrix
+    lambda1 = v_mean_normal * norm_
+    lambda2 = (v_mean_normal - c_mean_global) * norm_
+    lambda3 = (v_mean_normal + c_mean_global) * norm_
+
+    lambda1_p = positive_part(lambda1)
+    lambda2_p = positive_part(lambda2)
+    lambda3_p = positive_part(lambda3)
+    lambda2p3_half_p = 0.5 * (lambda2_p + lambda3_p)
+    lambda3m2_half_p = 0.5 * (lambda3_p - lambda2_p)
+
+    lambda1_m = negative_part(lambda1)
+    lambda2_m = negative_part(lambda2)
+    lambda3_m = negative_part(lambda3)
+    lambda2p3_half_m = 0.5 * (lambda2_m + lambda3_m)
+    lambda3m2_half_m = 0.5 * (lambda3_m - lambda2_m)
+
+    f1p = (lambda1_p * rho_prime_ll +
+           lambda3m2_half_p / c_mean_global * rho_mean_global * v_prime_normal_ll +
+           (lambda2p3_half_p - lambda1_p) / c_mean_global^2 * p_prime_ll)
+    f2p = (((lambda1_p * normal_vector[2]^2 +
+             lambda2p3_half_p * normal_vector[1]^2) * v1_prime_ll +
+            (lambda2p3_half_p - lambda1_p) * prod(normal_vector) * v2_prime_ll) +
+           lambda3m2_half_p / c_mean_global * normal_vector[1] * p_prime_ll /
+           rho_mean_global)
+    f3p = (((lambda1_p * normal_vector[1]^2 +
+             lambda2p3_half_p * normal_vector[2]^2) * v2_prime_ll +
+            (lambda2p3_half_p - lambda1_p) * prod(normal_vector) * v1_prime_ll) +
+           lambda3m2_half_p / c_mean_global * normal_vector[2] * p_prime_ll /
+           rho_mean_global)
+    f4p = (lambda3m2_half_p * c_mean_global * rho_mean_global * v_prime_normal_ll +
+           lambda2p3_half_p * p_prime_ll)
+
+    f1m = (lambda1_m * rho_prime_rr +
+           lambda3m2_half_m / c_mean_global * rho_mean_global * v_prime_normal_rr +
+           (lambda2p3_half_m - lambda1_m) / c_mean_global^2 * p_prime_rr)
+    f2m = (((lambda1_m * normal_vector[2]^2 +
+             lambda2p3_half_m * normal_vector[1]^2) * v1_prime_rr +
+            (lambda2p3_half_m - lambda1_m) * prod(normal_vector) * v2_prime_rr) +
+           lambda3m2_half_m / c_mean_global * normal_vector[1] * p_prime_rr /
+           rho_mean_global)
+    f3m = (((lambda1_m * normal_vector[1]^2 +
+             lambda2p3_half_m * normal_vector[2]^2) * v2_prime_rr +
+            (lambda2p3_half_m - lambda1_m) * prod(normal_vector) * v1_prime_rr) +
+           lambda3m2_half_m / c_mean_global * normal_vector[2] * p_prime_rr /
+           rho_mean_global)
+    f4m = (lambda3m2_half_m * c_mean_global * rho_mean_global * v_prime_normal_rr +
+           lambda2p3_half_m * p_prime_rr)
+
+    f1 = f1p + f1m
+    f2 = f2p + f2m
+    f3 = f3p + f3m
+    f4 = f4p + f4m
+
+    return SVector(f1, f2, f3, f4)
+end
+
 # Convert conservative variables to primitive
 @inline cons2prim(u, equations::LinearizedEulerEquations2D) = u
 @inline cons2entropy(u, ::LinearizedEulerEquations2D) = u