diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..05c09d57530 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,146 @@ +using OrdinaryDiffEq +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +cells_per_dimension = (64, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..f2e14273ae7 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,150 @@ +using OrdinaryDiffEq +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +# However TreeMesh will generate a 20_000 x 20_000 square domain instead +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index b0fd5c53f45..3c6f759db2b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -809,6 +809,98 @@ end return SVector(f1m, f2m, f3m, f4m) end +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" +struct FluxLMARS{SpeedOfSound} + # Estimate for the speed of sound + speed_of_sound::SpeedOfSound +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4 = v * u_ll + f4 = f4 + p_ll * v + else + f1, f2, f3, f4 = v * u_rr + f4 = f4 + p_rr * v + end + + if orientation == 1 + f2 = f2 + p + else # orientation == 2 + f3 = f3 + p + end + + return SVector(f1, f2, f3, f4) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4 = u_ll * v + f4 = f4 + p_ll * v + else + f1, f2, f3, f4 = u_rr * v + f4 = f4 + p_rr * v + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4) +end + """ splitting_vanleer_haenel(u, orientation::Integer, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 82c4a7efa32..292b912f009 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -944,11 +944,6 @@ References: Lagrangian Coordinate [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) """ -struct FluxLMARS{SpeedOfSound} - # Estimate for the speed of sound - speed_of_sound::SpeedOfSound -end - @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations3D) c = flux_lmars.speed_of_sound @@ -972,10 +967,14 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end if orientation == 1 @@ -985,7 +984,6 @@ end else # orientation == 3 f4 += p end - f5 += p * v return SVector(f1, f2, f3, f4, f5) end @@ -1011,18 +1009,21 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end - f2 += p * normal_direction[1] - f3 += p * normal_direction[2] - f4 += p * normal_direction[3] - f5 += p * v - - return SVector(f1, f2, f3, f4, f5) + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4 + p * normal_direction[3], + f5) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 4a2d2112c99..dc5d32b5a04 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -380,18 +380,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_circular_wind_nonconforming.jl"), l2=[ - 1.573832094977477e-7, - 3.863090659429634e-5, - 3.867293305754584e-5, - 3.686550296950078e-5, - 0.05508968493733932, + 1.5737711609657832e-7, + 3.8630261900166194e-5, + 3.8672287531936816e-5, + 3.6865116098660796e-5, + 0.05508620970403884, ], linf=[ - 2.2695202613887133e-6, - 0.0005314968179916946, - 0.0005314969614147458, - 0.0005130280733059617, - 0.7944959432352334, + 2.268845333053271e-6, + 0.000531462302113539, + 0.0005314624461298934, + 0.0005129931254772464, + 0.7942778058932163, ], tspan=(0.0, 2e2), coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI @@ -409,18 +409,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_baroclinic_instability.jl"), l2=[ - 6.725065410642336e-7, - 0.00021710117340245454, - 0.000438679759422352, - 0.00020836356588024185, - 0.07602006689579247, + 6.725093801700048e-7, + 0.00021710076010951073, + 0.0004386796338203878, + 0.00020836270267103122, + 0.07601887903440395, ], linf=[ - 1.9101671995258585e-5, - 0.029803626911022396, - 0.04847630924006063, - 0.022001371349740104, - 4.847761006938526, + 1.9107530539574924e-5, + 0.02980358831035801, + 0.048476331898047564, + 0.02200137344113612, + 4.848310144356219, ], tspan=(0.0, 1e2), # Decrease tolerance of adaptive time stepping to get similar results across different systems diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1addc29e3e6..e5d45ebcc07 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -611,6 +611,33 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_warm_bubble.jl"), + l2=[ + 0.00019387402388722496, + 0.03086514388623955, + 0.04541427917165, + 43.892826583444716, + ], + linf=[ + 0.0015942305974430138, + 0.17449778969139373, + 0.3729704262394843, + 307.6706958565337, + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.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)) < 100 + end +end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2=[ diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 65899cd5263..61b5c54b5e9 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -834,6 +834,32 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"), + l2=[ + 0.0001379946769624388, + 0.02078779689715382, + 0.033237241571263176, + 31.36068872331705, + ], + linf=[ + 0.0016286690573188434, + 0.15623770697198225, + 0.3341371832270615, + 334.5373488726036, + ], + tspan=(0.0, 10.0), + initial_refinement_level=4) + # 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)) < 100 + end +end + # Coverage test for all initial conditions @testset "Compressible Euler: Tests for initial conditions" begin @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin diff --git a/test/test_unit.jl b/test/test_unit.jl index 817b4cd550d..e8a8effbe29 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1287,6 +1287,49 @@ end end end +@timed_testset "Consistency check for LMARS flux" begin + equations = CompressibleEulerEquations2D(1.4) + flux_lmars = FluxLMARS(340) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + orientations = [1, 2] + u_values = [SVector(1.0, 0.5, -0.7, 1.0), + SVector(1.5, -0.2, 0.1, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + equations = CompressibleEulerEquations3D(1.4) + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + orientations = [1, 2, 3] + u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), + SVector(1.5, -0.2, 0.1, 0.2, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + @testset "FluxRotated vs. direct implementation" begin @timed_testset "CompressibleEulerMulticomponentEquations2D" begin equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), @@ -1320,7 +1363,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std)