1D Linearized Euler

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

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

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,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()

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