diff --git a/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl b/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl new file mode 100644 index 0000000000..fd066ef895 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl @@ -0,0 +1,64 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations3D(v_mean_global = (0.0, 0.0, 0.0), c_mean_global = 1.0, + rho_mean_global = 1.0) + +initial_condition = initial_condition_convergence_test + +solver = DGSEM(polydeg = 3, surface_flux = flux_hll) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# `initial_refinement_level` is provided here to allow for a +# convenient convergence test, see +# https://trixi-framework.github.io/Trixi.jl/stable/#Performing-a-convergence-analysis +trees_per_dimension = (4, 4, 4) +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 0.2 +tspan = (0.0, 0.2) +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() + +analysis_interval = 100 + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.8) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + 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() # print the timer summary diff --git a/examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl new file mode 100644 index 0000000000..2eb8ef822a --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -0,0 +1,65 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations3D(v_mean_global = (0.5, 0.5, 0.5), c_mean_global = 1.0, + rho_mean_global = 1.0) + +solver = DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (90.0, 90.0, 90.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000, + periodicity = false) + +# Initialize density and pressure perturbation with a Gaussian bump +# that splits into radial waves which are advected with v - c and v + c. +function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations3D) + v1_prime = 0.0 + v2_prime = 0.0 + v3_prime = 0.0 + rho_prime = p_prime = 2 * exp(-((x[1] - 45)^2 + (x[2] - 45)^2) / 25) + return SVector(rho_prime, v1_prime, v2_prime, v3_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. + +# At t = 30, the wave moving with v + c crashes into the wall +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 StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.9) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, + 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/Trixi.jl b/src/Trixi.jl index 3e7fc5ee51..5511f7e02e 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -160,7 +160,7 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations1D, LinearizedEulerEquations2D, + LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D, PolytropicEulerEquations2D, TrafficFlowLWREquations1D diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ef4dac8b1a..69a54f7faf 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -503,6 +503,7 @@ abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("linearized_euler_1d.jl") include("linearized_euler_2d.jl") +include("linearized_euler_3d.jl") abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <: AbstractEquations{NDIMS, NVARS} end diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 856a53a3d5..19a3bdcb3b 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -8,7 +8,7 @@ @doc raw""" LinearizedEulerEquations1D(v_mean_global, c_mean_global, rho_mean_global) -Linearized euler equations in one space dimension. The equations are given by +Linearized Euler equations in one space dimension. The equations are given by ```math \partial_t \begin{pmatrix} diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index d9d0d44e55..3df3093069 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/src/equations/linearized_euler_2d.jl @@ -8,7 +8,7 @@ @doc raw""" LinearizedEulerEquations2D(v_mean_global, c_mean_global, rho_mean_global) -Linearized euler equations in two space dimensions. The equations are given by +Linearized Euler equations in two space dimensions. The equations are given by ```math \partial_t \begin{pmatrix} diff --git a/src/equations/linearized_euler_3d.jl b/src/equations/linearized_euler_3d.jl new file mode 100644 index 0000000000..ab5f9863db --- /dev/null +++ b/src/equations/linearized_euler_3d.jl @@ -0,0 +1,254 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + LinearizedEulerEquations3D(v_mean_global, c_mean_global, rho_mean_global) + +Linearized Euler equations in three space dimensions. The equations are given by +```math +\partial_t +\begin{pmatrix} + \rho' \\ v_1' \\ v_2' \\ v_3' \\ p' +\end{pmatrix} ++ +\partial_x +\begin{pmatrix} + \bar{\rho} v_1' + \bar{v_1} \rho ' \\ + \bar{v_1} v_1' + \frac{p'}{\bar{\rho}} \\ + \bar{v_1} v_2' \\ + \bar{v_1} v_3' \\ + \bar{v_1} p' + c^2 \bar{\rho} v_1' +\end{pmatrix} ++ +\partial_y +\begin{pmatrix} + \bar{\rho} v_2' + \bar{v_2} \rho ' \\ + \bar{v_2} v_1' \\ + \bar{v_2} v_2' + \frac{p'}{\bar{\rho}} \\ + \bar{v_2} v_3' \\ + \bar{v_2} p' + c^2 \bar{\rho} v_2' +\end{pmatrix} ++ +\partial_z +\begin{pmatrix} + \bar{\rho} v_3' + \bar{v_3} \rho ' \\ + \bar{v_3} v_1' \\ + \bar{v_3} v_2' \\ + \bar{v_3} v_3' + \frac{p'}{\bar{\rho}} \\ + \bar{v_3} p' + c^2 \bar{\rho} v_3' +\end{pmatrix} += +\begin{pmatrix} + 0 \\ 0 \\ 0 \\ 0 \\ 0 +\end{pmatrix} +``` +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and ``c`` is the speed of sound. +The unknowns are the acoustic velocities ``v' = (v_1', v_2, v_3')``, the pressure ``p'`` and the density ``\rho'``. +""" +struct LinearizedEulerEquations3D{RealT <: Real} <: + AbstractLinearizedEulerEquations{3, 5} + v_mean_global::SVector{3, RealT} + c_mean_global::RealT + rho_mean_global::RealT +end + +function LinearizedEulerEquations3D(v_mean_global::NTuple{3, <:Real}, + c_mean_global::Real, rho_mean_global::Real) + if rho_mean_global < 0 + throw(ArgumentError("rho_mean_global must be non-negative")) + elseif c_mean_global < 0 + throw(ArgumentError("c_mean_global must be non-negative")) + end + + return LinearizedEulerEquations3D(SVector(v_mean_global), c_mean_global, + rho_mean_global) +end + +function LinearizedEulerEquations3D(; v_mean_global::NTuple{3, <:Real}, + c_mean_global::Real, rho_mean_global::Real) + return LinearizedEulerEquations3D(v_mean_global, c_mean_global, + rho_mean_global) +end + +function varnames(::typeof(cons2cons), ::LinearizedEulerEquations3D) + ("rho_prime", "v1_prime", "v2_prime", "v3_prime", "p_prime") +end +function varnames(::typeof(cons2prim), ::LinearizedEulerEquations3D) + ("rho_prime", "v1_prime", "v2_prime", "v3_prime", "p_prime") +end + +""" + initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations3D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations3D) + rho_prime = -cospi(2 * t) * (sinpi(2 * x[1]) + sinpi(2 * x[2]) + sinpi(2 * x[3])) + v1_prime = sinpi(2 * t) * cospi(2 * x[1]) + v2_prime = sinpi(2 * t) * cospi(2 * x[2]) + v3_prime = sinpi(2 * t) * cospi(2 * x[3]) + p_prime = rho_prime + + return SVector(rho_prime, v1_prime, v2_prime, v3_prime, p_prime) +end + +""" + boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function, + equations::LinearizedEulerEquations3D) + +Boundary conditions for a solid wall. +""" +function boundary_condition_wall(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::LinearizedEulerEquations3D) + # Boundary state is equal to the inner state except for the velocity. For boundaries + # in the -x/+x direction, we multiply the velocity in the x direction by -1. + # Similarly, for boundaries in the -y/+y or -z/+z direction, we multiply the + # velocity in the y or z direction by -1 + if direction in (1, 2) # x direction + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4], + u_inner[5]) + elseif direction in (3, 4) # y direction + u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4], + u_inner[5]) + else # z direction = (5, 6) + u_boundary = SVector(u_inner[1], u_inner[2], u_inner[3], -u_inner[4], + u_inner[5]) + end + + # Calculate boundary flux depending on the orientation of the boundary + # Odd directions are in negative coordinate direction, + # even directions are in positive coordinate direction. + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# Calculate 3D flux for a single point +@inline function flux(u, orientation::Integer, equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, v2_prime, v3_prime, p_prime = u + if orientation == 1 + f1 = v_mean_global[1] * rho_prime + rho_mean_global * v1_prime + f2 = v_mean_global[1] * v1_prime + p_prime / rho_mean_global + f3 = v_mean_global[1] * v2_prime + f4 = v_mean_global[1] * v3_prime + f5 = v_mean_global[1] * p_prime + c_mean_global^2 * rho_mean_global * v1_prime + elseif orientation == 2 + f1 = v_mean_global[2] * rho_prime + rho_mean_global * v2_prime + f2 = v_mean_global[2] * v1_prime + f3 = v_mean_global[2] * v2_prime + p_prime / rho_mean_global + f4 = v_mean_global[2] * v3_prime + f5 = v_mean_global[2] * p_prime + c_mean_global^2 * rho_mean_global * v2_prime + else # orientation == 3 + f1 = v_mean_global[3] * rho_prime + rho_mean_global * v3_prime + f2 = v_mean_global[3] * v1_prime + f3 = v_mean_global[3] * v2_prime + f4 = v_mean_global[3] * v3_prime + p_prime / rho_mean_global + f5 = v_mean_global[3] * p_prime + c_mean_global^2 * rho_mean_global * v3_prime + end + + return SVector(f1, f2, f3, f4, f5) +end + +# Calculate 3D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, v2_prime, v3_prime, p_prime = u + + v_mean_normal = v_mean_global[1] * normal_direction[1] + + v_mean_global[2] * normal_direction[2] + + v_mean_global[3] * normal_direction[3] + v_prime_normal = v1_prime * normal_direction[1] + v2_prime * normal_direction[2] + + v3_prime * normal_direction[3] + + 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 * v3_prime + normal_direction[3] * p_prime / rho_mean_global + f5 = v_mean_normal * p_prime + c_mean_global^2 * rho_mean_global * v_prime_normal + + return SVector(f1, f2, f3, f4, f5) +end + +@inline have_constant_speed(::LinearizedEulerEquations3D) = True() + +@inline function max_abs_speeds(equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global[1]) + c_mean_global, abs(v_mean_global[2]) + c_mean_global, + abs(v_mean_global[3]) + c_mean_global +end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + if orientation == 1 + return abs(v_mean_global[1]) + c_mean_global + elseif orientation == 2 + return abs(v_mean_global[2]) + c_mean_global + else # orientation == 3 + return abs(v_mean_global[3]) + c_mean_global + end +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @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] + + normal_direction[3] * v_mean_global[3] + return abs(v_mean_normal) + c_mean_global * norm(normal_direction) +end + +# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations3D) + min_max_speed_davis(u_ll, u_rr, orientation, equations) +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + min_max_speed_davis(u_ll, u_rr, normal_direction, equations) +end + +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + + λ_min = v_mean_global[orientation] - c_mean_global + λ_max = v_mean_global[orientation] + c_mean_global + + return λ_min, λ_max +end + +@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + + norm_ = norm(normal_direction) + + v_normal = v_mean_global[1] * normal_direction[1] + + v_mean_global[2] * normal_direction[2] + + v_mean_global[3] * normal_direction[3] + + # The v_normals are already scaled by the norm + λ_min = v_normal - c_mean_global * norm_ + λ_max = v_normal + c_mean_global * norm_ + + return λ_min, λ_max +end + +# Convert conservative variables to primitive +@inline cons2prim(u, equations::LinearizedEulerEquations3D) = u +@inline cons2entropy(u, ::LinearizedEulerEquations3D) = u +end # muladd diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index ea7d9193ad..8b684b22d0 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -535,6 +535,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_linearizedeuler_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_linearizedeuler_convergence.jl"), + l2=[ + 0.04452389418193219, 0.03688186699434862, + 0.03688186699434861, 0.03688186699434858, + 0.044523894181932186, + ], + linf=[ + 0.2295447498696467, 0.058369658071546704, + 0.05836965807154648, 0.05836965807154648, 0.2295447498696468, + ]) + # 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 end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_tree_3d_linearizedeuler.jl b/test/test_tree_3d_linearizedeuler.jl new file mode 100644 index 0000000000..00f8d62dad --- /dev/null +++ b/test/test_tree_3d_linearizedeuler.jl @@ -0,0 +1,35 @@ + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem") + +@testset "Linearized Euler Equations 3D" begin +#! format: noindent + +@trixi_testset "elixir_linearizedeuler_gauss_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gauss_wall.jl"), + l2=[ + 0.020380328336745232, 0.027122442311921492, + 0.02712244231192152, 8.273108096127844e-17, + 0.020380328336745232, + ], + linf=[ + 0.2916021983572774, 0.32763703462270843, + 0.32763703462270855, 1.641012595221666e-15, + 0.2916021983572774, + ], + tspan=(0.0, 1.0)) + + # 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 +end diff --git a/test/test_tree_3d_part2.jl b/test/test_tree_3d_part2.jl index 5c7301654a..4b9da039f9 100644 --- a/test/test_tree_3d_part2.jl +++ b/test/test_tree_3d_part2.jl @@ -22,6 +22,9 @@ isdir(outdir) && rm(outdir, recursive = true) # Compressible Euler with self-gravity include("test_tree_3d_eulergravity.jl") + + # Linearized Euler + include("test_tree_3d_linearizedeuler.jl") end @trixi_testset "Additional tests in 3D" begin diff --git a/test/test_unit.jl b/test/test_unit.jl index 9a4ae2ede1..b23d4aa60e 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1536,6 +1536,60 @@ end end end +@testset "Equivalent Wave Speed Estimates" begin + @timed_testset "Linearized Euler 3D" begin + equations = LinearizedEulerEquations3D(v_mean_global = (0.42, 0.37, 0.7), + c_mean_global = 1.0, + rho_mean_global = 1.0) + + normal_x = SVector(1.0, 0.0, 0.0) + normal_y = SVector(0.0, 1.0, 0.0) + normal_z = SVector(0.0, 0.0, 1.0) + + u_ll = SVector(0.3, 0.5, -0.7, 0.1, 1.0) + u_rr = SVector(0.5, -0.2, 0.1, 0.2, 5.0) + + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 1, equations), + max_abs_speed_naive(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 2, equations), + max_abs_speed_naive(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 3, equations), + max_abs_speed_naive(u_ll, u_rr, normal_z, + equations))) + + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 1, equations), + min_max_speed_naive(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 2, equations), + min_max_speed_naive(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 3, equations), + min_max_speed_naive(u_ll, u_rr, normal_z, + equations))) + + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 1, equations), + min_max_speed_davis(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 2, equations), + min_max_speed_davis(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 3, equations), + min_max_speed_davis(u_ll, u_rr, normal_z, + equations))) + end +end + @testset "SimpleKronecker" begin N = 3