From 5ff677c1d246e7a500ab845201bc328d1d3bde92 Mon Sep 17 00:00:00 2001 From: Lars Christmann Date: Tue, 11 Jul 2023 18:53:25 +0200 Subject: [PATCH] Implement upwind flux for linearized Euler equations (#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 --------- Co-authored-by: Michael Schlottke-Lakemper --- .../elixir_linearizedeuler_gaussian_source.jl | 89 ++++++++ .../elixir_linearizedeuler_gauss_wall.jl | 68 ++++++ src/equations/linearized_euler_2d.jl | 212 +++++++++++++++++- test/test_p4est_2d.jl | 5 + test/test_tree_2d_linearizedeuler.jl | 6 + test/test_unit.jl | 20 ++ 6 files changed, 399 insertions(+), 1 deletion(-) create mode 100644 examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl create mode 100644 examples/tree_2d_dgsem/elixir_linearizedeuler_gauss_wall.jl diff --git a/examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl b/examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl new file mode 100644 index 0000000000..ba2ec82777 --- /dev/null +++ b/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() diff --git a/examples/tree_2d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_2d_dgsem/elixir_linearizedeuler_gauss_wall.jl new file mode 100644 index 0000000000..14fe201a29 --- /dev/null +++ b/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() diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index cd681365ca..e478c32bd2 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/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 diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f66664c7a8..c4ce2619e1 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -164,6 +164,11 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.02)) end + @trixi_testset "elixir_linearizedeuler_gaussian_source.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gaussian_source.jl"), + l2 = [0.006047938590548741, 0.0040953286019907035, 0.004222698522497298, 0.006269492499336128], + linf = [0.06386175207349379, 0.0378926444850457, 0.041759728067967065, 0.06430136016259067]) + end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_tree_2d_linearizedeuler.jl b/test/test_tree_2d_linearizedeuler.jl index 540b395121..2c5f6dc2cd 100644 --- a/test/test_tree_2d_linearizedeuler.jl +++ b/test/test_tree_2d_linearizedeuler.jl @@ -13,4 +13,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [0.0011006084408365924, 0.0005788678074691855, 0.0005788678074701847, 0.0011006084408365924] ) end + + @trixi_testset "elixir_linearizedeuler_gauss_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gauss_wall.jl"), + l2 = [0.048185623945503485, 0.01941899333212175, 0.019510224816991825, 0.048185623945503485], + linf = [1.0392165942153189, 0.18188777290819994, 0.1877028372108587, 1.0392165942153189]) + end end diff --git a/test/test_unit.jl b/test/test_unit.jl index 2156e9bac3..b0c3e4205e 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -670,6 +670,26 @@ isdir(outdir) && rm(outdir, recursive=true) for normal_direction in normal_directions @test flux_godunov(u, u, normal_direction, equation) ≈ flux(u, normal_direction, equation) end + + # Linearized Euler 2D + equation = LinearizedEulerEquations2D(v_mean_global=(0.5, -0.7), c_mean_global=1.1, + rho_mean_global=1.2) + u_values = [SVector(1.0, 0.5, -0.7, 1.0), + SVector(1.5, -0.2, 0.1, 5.0),] + + orientations = [1, 2] + for orientation in orientations, u in u_values + @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions, u in u_values + @test flux_godunov(u, u, normal_direction, equation) ≈ flux(u, normal_direction, equation) + end end @timed_testset "Consistency check for Engquist-Osher flux" begin