From 5b6eead5006bf2790b5afc776165d589ef066544 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 8 Mar 2024 15:51:04 +0100 Subject: [PATCH 01/12] 1D Linearized Euler --- ...r_linearizedeuler_characteristic_system.jl | 109 +++++++++++++ .../elixir_linearizedeuler_convergence.jl | 60 ++++++++ .../elixir_linearizedeuler_gauss_wall.jl | 63 ++++++++ src/Trixi.jl | 2 +- src/equations/equations.jl | 1 + src/equations/linearized_euler_1d.jl | 143 ++++++++++++++++++ test/test_structured_1d.jl | 15 ++ test/test_tree_1d.jl | 3 + test/test_tree_1d_linearizedeuler.jl | 52 +++++++ 9 files changed, 447 insertions(+), 1 deletion(-) create mode 100644 examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl create mode 100644 examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl create mode 100644 examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl create mode 100644 src/equations/linearized_euler_1d.jl create mode 100644 test/test_tree_1d_linearizedeuler.jl 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..e26d10446f3 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using LinearAlgebra: dot +using Trixi + +############################################################################### +# semidiscretization of the 1 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 +LinEuler_EigVals = [v_0 - c_0; v_0; v_0 + c_0] +LinEuler_EigVecs = [-rho_0/c_0 1 rho_0/c_0; + 1 0 1; + -rho_0*c_0 0 rho_0*c_0] +LinEuler_EigVecs_inv = inv(LinEuler_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 * LinEuler_EigVals +end + +function compute_primal_sol(char_vars) + return LinEuler_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) + for p in 1:3 + w[p] = dot(LinEuler_EigVecs_inv[p, :], + initial_condition_entropy_wave(x_char[p], 0, equations)) # Assumes t_0 = 0 + 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..7fbe3e4d231 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl @@ -0,0 +1,60 @@ +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); + +# print the timer summary +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..c5f06fd0259 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -0,0 +1,63 @@ + +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,) + +# Create a uniformly refined mesh with periodic boundaries +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() diff --git a/src/Trixi.jl b/src/Trixi.jl index da7359999c5..1baa5eae806 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -161,7 +161,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..c5e539d7028 --- /dev/null +++ b/src/equations/linearized_euler_1d.jl @@ -0,0 +1,143 @@ +# 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 \\ 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')``, 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 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..ad8181bf5ce --- /dev/null +++ b/test/test_tree_1d_linearizedeuler.jl @@ -0,0 +1,52 @@ + +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 From 21536b2dc12f9dbca4f69856b73911eb33b0804b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 8 Mar 2024 15:56:43 +0100 Subject: [PATCH 02/12] remove redundant comment --- examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl index 7fbe3e4d231..5b17ab4f3dc 100644 --- a/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl @@ -56,5 +56,4 @@ 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 From 0a5398afd3a1ec21d642949780197477a30e2c1d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 9 Mar 2024 16:54:34 +0100 Subject: [PATCH 03/12] Fix docu --- src/equations/linearized_euler_1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index c5e539d7028..27245a4ddd3 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -12,7 +12,7 @@ Linearized euler equations in one space dimension. The equations are given by ```math \partial_t \begin{pmatrix} - \rho' \\ v_1' \\ \\ p' + \rho' \\ v_1' \\ p' \end{pmatrix} + \partial_x @@ -21,11 +21,11 @@ Linearized euler equations in one space dimension. The equations are given by \end{pmatrix} = \begin{pmatrix} - 0 \\ 0 \\ 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')``, the pressure ``p'`` and the density ``\rho'``. +The unknowns are the acoustic velocity ``v_1'``, the pressure ``p'`` and the density ``\rho'``. """ struct LinearizedEulerEquations1D{RealT <: Real} <: AbstractLinearizedEulerEquations{1, 3} From 651c63d35dadb5f7592871434f30e19aad3d8119 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 11 Mar 2024 09:57:03 +0100 Subject: [PATCH 04/12] Update src/equations/linearized_euler_1d.jl --- src/equations/linearized_euler_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 27245a4ddd3..559fb1dc0f0 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -46,7 +46,7 @@ function LinearizedEulerEquations1D(v_mean_global::Real, rho_mean_global) end -# Constructor with keywords (note the leading ';) +# 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, From ebe0a7b9b62d18628df6e936ff952ed4b9792d22 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 12 Mar 2024 15:15:39 +0100 Subject: [PATCH 05/12] Update examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl --- .../elixir_linearizedeuler_characteristic_system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl index e26d10446f3..c1a11f990a8 100644 --- a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -4,7 +4,7 @@ using LinearAlgebra: dot using Trixi ############################################################################### -# semidiscretization of the 1 linearized Euler equations +# semidiscretization of the linearized Euler equations rho_0 = 1.0 v_0 = 1.0 From 287e15bc044f3bd54d9c529037acb6f9794b1dd1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 12 Mar 2024 15:27:01 +0100 Subject: [PATCH 06/12] code review --- ...lixir_linearizedeuler_characteristic_system.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl index e26d10446f3..068df61e00d 100644 --- a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -20,20 +20,20 @@ cells_per_dimension = (64,) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) # Linearized Euler: Eigensystem -LinEuler_EigVals = [v_0 - c_0; v_0; v_0 + c_0] -LinEuler_EigVecs = [-rho_0/c_0 1 rho_0/c_0; +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] -LinEuler_EigVecs_inv = inv(LinEuler_EigVecs) +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 * LinEuler_EigVals + return SVector(x[1], x[1], x[1]) .- t * lin_euler_eigvals end function compute_primal_sol(char_vars) - return LinEuler_EigVecs * char_vars + return lin_euler_eigvecs * char_vars end # Initial condition is in principle arbitrary, only periodicity is required @@ -66,9 +66,10 @@ function initial_condition_char_vars(x, t, equations::LinearizedEulerEquations1D # Set up characteristic variables w = zeros(3) + t_0 = 0 # Assumes t_0 = 0 for p in 1:3 - w[p] = dot(LinEuler_EigVecs_inv[p, :], - initial_condition_entropy_wave(x_char[p], 0, equations)) # Assumes t_0 = 0 + 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) From 7cff25b15bee21e9ec73f7ec99083b8b210ca590 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 12 Mar 2024 15:36:15 +0100 Subject: [PATCH 07/12] fmt --- .../elixir_linearizedeuler_characteristic_system.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl index e126134a655..29f4bd0a8cf 100644 --- a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -22,8 +22,8 @@ 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] + 1 0 1; + -rho_0*c_0 0 rho_0*c_0] lin_euler_eigvecs_inv = inv(lin_euler_eigvecs) # Trace back characteristics. @@ -68,7 +68,7 @@ function initial_condition_char_vars(x, t, equations::LinearizedEulerEquations1D 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) + u_char = initial_condition_entropy_wave(x_char[p], t_0, equations) w[p] = dot(lin_euler_eigvecs_inv[p, :], u_char) end From c71fe7ddff28b5772854941c36c526cdb2742e3d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 28 Mar 2024 08:23:47 +0100 Subject: [PATCH 08/12] Apply suggestions from code review Co-authored-by: Andrew Winters --- examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl | 1 - src/equations/linearized_euler_1d.jl | 3 ++- test/test_tree_1d_linearizedeuler.jl | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl index c5f06fd0259..c087989e534 100644 --- a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -13,7 +13,6 @@ solver = DGSEM(polydeg = 5, surface_flux = flux_hll) coordinates_min = (0.0,) coordinates_max = (90.0,) -# Create a uniformly refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 6, n_cells_max = 100_000, diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 559fb1dc0f0..04ea812a554 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -25,7 +25,8 @@ Linearized euler equations in one space dimension. The equations are given by \end{pmatrix} ``` The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. -The unknowns are the acoustic velocity ``v_1'``, the pressure ``p'`` and the density ``\rho'``. +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} diff --git a/test/test_tree_1d_linearizedeuler.jl b/test/test_tree_1d_linearizedeuler.jl index ad8181bf5ce..c7cffee3f66 100644 --- a/test/test_tree_1d_linearizedeuler.jl +++ b/test/test_tree_1d_linearizedeuler.jl @@ -39,7 +39,6 @@ end 0.9999720404625275, 1.9999505145390108, ]) - # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From aa2d8e34a55230281b792a9c9ea29af530b1d412 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 28 Mar 2024 08:28:00 +0100 Subject: [PATCH 09/12] Update src/equations/linearized_euler_1d.jl --- src/equations/linearized_euler_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 04ea812a554..22f452219b5 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -47,7 +47,7 @@ function LinearizedEulerEquations1D(v_mean_global::Real, rho_mean_global) end -# Constructor with keywords (note the leading ';') +# 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, From ec0f3b05fde492da99d6efe9e5eb71a1b85bee73 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 28 Mar 2024 15:16:02 +0100 Subject: [PATCH 10/12] Comments --- examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl index c087989e534..0884249559a 100644 --- a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -18,6 +18,9 @@ mesh = TreeMesh(coordinates_min, coordinates_max, 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) From 5d1e2b2cb436b1defec15037d22489e76d7b0f13 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 28 Mar 2024 15:28:07 +0100 Subject: [PATCH 11/12] comments --- .../elixir_linearizedeuler_characteristic_system.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl index 29f4bd0a8cf..663b25b18c0 100644 --- a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -19,6 +19,8 @@ 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; From ea1b268b952d1fcf722017ad673191c00f63456b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 28 Mar 2024 15:29:19 +0100 Subject: [PATCH 12/12] news --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 022252e61a9..148bd5c246f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ for human readability. #### Added - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. - +- Implementation of 1D Linearized Euler Equations. ## Changes when updating to v0.7 from v0.6.x