1D Linearized Euler

examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl

@@ -0,0 +1,110 @@
+
+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)
+
+# 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

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

examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl

@@ -0,0 +1,62 @@
+
+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)
+
+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()

src/Trixi.jl

@@ -160,7 +160,7 @@ export AcousticPerturbationEquations2D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
        ShallowWaterEquationsQuasi1D,
-       LinearizedEulerEquations2D,
+       LinearizedEulerEquations1D, LinearizedEulerEquations2D,
        PolytropicEulerEquations2D,
        TrafficFlowLWREquations1D
 

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} <:

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 (note the leading ';')
+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