diff --git a/NEWS.md b/NEWS.md index 5516fc9add5..e1238562fc9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,7 @@ for human readability. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. +- Implementation of 1D Linearized Euler Equations. - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. ## Changes when updating to v0.7 from v0.6.x diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl new file mode 100644 index 00000000000..663b25b18c0 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -0,0 +1,112 @@ + +using OrdinaryDiffEq +using LinearAlgebra: dot +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +rho_0 = 1.0 +v_0 = 1.0 +c_0 = 1.0 +equations = LinearizedEulerEquations1D(rho_0, v_0, c_0) + +solver = DGSEM(polydeg = 3, surface_flux = flux_hll) + +coordinates_min = (0.0,) # minimum coordinate +coordinates_max = (1.0,) # maximum coordinate +cells_per_dimension = (64,) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# For eigensystem of the linearized Euler equations see e.g. +# https://www.nas.nasa.gov/assets/nas/pdf/ams/2018/introtocfd/Intro2CFD_Lecture1_Pulliam_Euler_WaveEQ.pdf +# Linearized Euler: Eigensystem +lin_euler_eigvals = [v_0 - c_0; v_0; v_0 + c_0] +lin_euler_eigvecs = [-rho_0/c_0 1 rho_0/c_0; + 1 0 1; + -rho_0*c_0 0 rho_0*c_0] +lin_euler_eigvecs_inv = inv(lin_euler_eigvecs) + +# Trace back characteristics. +# See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf, p.95 +function compute_char_initial_pos(x, t) + return SVector(x[1], x[1], x[1]) .- t * lin_euler_eigvals +end + +function compute_primal_sol(char_vars) + return lin_euler_eigvecs * char_vars +end + +# Initial condition is in principle arbitrary, only periodicity is required +function initial_condition_entropy_wave(x, t, equations::LinearizedEulerEquations1D) + # Parameters + alpha = 1.0 + beta = 150.0 + center = 0.5 + + rho_prime = alpha * exp(-beta * (x[1] - center)^2) + v_prime = 0.0 + p_prime = 0.0 + + return SVector(rho_prime, v_prime, p_prime) +end + +function initial_condition_char_vars(x, t, equations::LinearizedEulerEquations1D) + # Trace back characteristics + x_char = compute_char_initial_pos(x, t) + + # Employ periodicity + for p in 1:3 + while x_char[p] < coordinates_min[1] + x_char[p] += coordinates_max[1] - coordinates_min[1] + end + while x_char[p] > coordinates_max[1] + x_char[p] -= coordinates_max[1] - coordinates_min[1] + end + end + + # Set up characteristic variables + w = zeros(3) + t_0 = 0 # Assumes t_0 = 0 + for p in 1:3 + u_char = initial_condition_entropy_wave(x_char[p], t_0, equations) + w[p] = dot(lin_euler_eigvecs_inv[p, :], u_char) + end + + return compute_primal_sol(w) +end + +initial_condition = initial_condition_char_vars + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.3) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +# 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 + +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); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl new file mode 100644 index 00000000000..5b17ab4f3dc --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl @@ -0,0 +1,59 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations1D(v_mean_global = 0.0, c_mean_global = 1.0, + rho_mean_global = 1.0) + +initial_condition = initial_condition_convergence_test + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-1.0) +coordinates_max = (1.0) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.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() + +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); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl new file mode 100644 index 00000000000..0884249559a --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -0,0 +1,65 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations1D(v_mean_global = 0.5, c_mean_global = 1.0, + rho_mean_global = 1.0) + +solver = DGSEM(polydeg = 5, surface_flux = flux_hll) + +coordinates_min = (0.0,) +coordinates_max = (90.0,) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 100_000, + periodicity = false) + +# Initialize density and pressure perturbation with a Gaussian bump +# that is advected to left with v - c and to the right with v + c. +# Correspondigly, the bump splits in half. +function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D) + v1_prime = 0.0 + rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25) + return SVector(rho_prime, v1_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 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, + 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 168441513ad..883f8d66f07 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -160,7 +160,7 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations2D, + LinearizedEulerEquations1D, LinearizedEulerEquations2D, PolytropicEulerEquations2D, TrafficFlowLWREquations1D diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8f476cf6f16..ef4dac8b1a5 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -501,6 +501,7 @@ include("acoustic_perturbation_2d.jl") # Linearized Euler equations abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("linearized_euler_1d.jl") include("linearized_euler_2d.jl") abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <: diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl new file mode 100644 index 00000000000..22f452219b5 --- /dev/null +++ b/src/equations/linearized_euler_1d.jl @@ -0,0 +1,144 @@ +# 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""" + LinearizedEulerEquations1D(v_mean_global, c_mean_global, rho_mean_global) + +Linearized euler equations in one space dimension. The equations are given by +```math +\partial_t +\begin{pmatrix} + \rho' \\ v_1' \\ 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} p' + c^2 \bar{\rho} v_1' +\end{pmatrix} += +\begin{pmatrix} + 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 perturbation quantities of the acoustic velocity ``v_1'``, the pressure ``p'`` +and the density ``\rho'``. +""" +struct LinearizedEulerEquations1D{RealT <: Real} <: + AbstractLinearizedEulerEquations{1, 3} + v_mean_global::RealT + c_mean_global::RealT + rho_mean_global::RealT +end + +function LinearizedEulerEquations1D(v_mean_global::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 LinearizedEulerEquations1D(v_mean_global, c_mean_global, + rho_mean_global) +end + +# Constructor with keywords +function LinearizedEulerEquations1D(; v_mean_global::Real, + c_mean_global::Real, rho_mean_global::Real) + return LinearizedEulerEquations1D(v_mean_global, c_mean_global, + rho_mean_global) +end + +function varnames(::typeof(cons2cons), ::LinearizedEulerEquations1D) + ("rho_prime", "v1_prime", "p_prime") +end +function varnames(::typeof(cons2prim), ::LinearizedEulerEquations1D) + ("rho_prime", "v1_prime", "p_prime") +end + +""" + initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations1D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations1D) + rho_prime = -cospi(2 * t) * sinpi(2 * x[1]) + v1_prime = sinpi(2 * t) * cospi(2 * x[1]) + p_prime = rho_prime + + return SVector(rho_prime, v1_prime, p_prime) +end + +""" + boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function, + equations::LinearizedEulerEquations1D) + +Boundary conditions for a solid wall. +""" +function boundary_condition_wall(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::LinearizedEulerEquations1D) + # 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. + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3]) + + # Calculate boundary flux + 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 1D flux for a single point +@inline function flux(u, orientation::Integer, equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, p_prime = u + f1 = v_mean_global * rho_prime + rho_mean_global * v1_prime + f2 = v_mean_global * v1_prime + p_prime / rho_mean_global + f3 = v_mean_global * p_prime + c_mean_global^2 * rho_mean_global * v1_prime + + return SVector(f1, f2, f3) +end + +@inline have_constant_speed(::LinearizedEulerEquations1D) = True() + +@inline function max_abs_speeds(equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global) + c_mean_global +end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global) + c_mean_global +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::LinearizedEulerEquations1D) + min_max_speed_davis(u_ll, u_rr, orientation, 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::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + + λ_min = v_mean_global - c_mean_global + λ_max = v_mean_global + c_mean_global + + return λ_min, λ_max +end + +# Convert conservative variables to primitive +@inline cons2prim(u, equations::LinearizedEulerEquations1D) = u +@inline cons2entropy(u, ::LinearizedEulerEquations1D) = u +end # muladd diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index fea06554c57..f97696d089a 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -139,6 +139,21 @@ end end end +@trixi_testset "elixir_linearizedeuler_characteristic_system.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_linearizedeuler_characteristic_system.jl"), + l2=[2.9318078842789714e-6, 0.0, 0.0], + linf=[4.291208715723194e-5, 0.0, 0.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 + @trixi_testset "elixir_traffic_flow_lwr_greenlight.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_traffic_flow_lwr_greenlight.jl"), diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 4a25a51a45e..76129f15e07 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -48,6 +48,9 @@ isdir(outdir) && rm(outdir, recursive = true) # Traffic flow LWR include("test_tree_1d_traffic_flow_lwr.jl") + + # Linearized Euler + include("test_tree_1d_linearizedeuler.jl") end # Coverage test for all initial conditions diff --git a/test/test_tree_1d_linearizedeuler.jl b/test/test_tree_1d_linearizedeuler.jl new file mode 100644 index 00000000000..c7cffee3f66 --- /dev/null +++ b/test/test_tree_1d_linearizedeuler.jl @@ -0,0 +1,51 @@ + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") + +@testset "Linearized Euler Equations 1D" begin +#! format: noindent + +@trixi_testset "elixir_linearizedeuler_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_convergence.jl"), + l2=[ + 0.00010894927270421941, + 0.00014295255695912358, + 0.00010894927270421941, + ], + linf=[ + 0.0005154647164193893, + 0.00048457837684242266, + 0.0005154647164193893, + ]) + # 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_linearizedeuler_gauss_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gauss_wall.jl"), + l2=[0.650082087850354, 0.2913911415488769, 0.650082087850354], + linf=[ + 1.9999505145390108, + 0.9999720404625275, + 1.9999505145390108, + ]) + # 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