From 641fde8c40e5b182193320fb805def4266d894d8 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 2 Apr 2024 09:31:47 +0200 Subject: [PATCH 01/17] Bump crate-ci/typos from 1.18.2 to 1.19.0 (#1894) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.18.2 to 1.19.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.18.2...v1.19.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 87e34cb50f3..e8f2271f831 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.18.2 + uses: crate-ci/typos@v1.19.0 From 09600b56915296e28e5792f87c91a820e9c4eee4 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 3 Apr 2024 08:15:22 +0200 Subject: [PATCH 02/17] Bump julia-actions/setup-julia from 1 to 2 (#1895) Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 1 to 2. - [Release notes](https://github.com/julia-actions/setup-julia/releases) - [Commits](https://github.com/julia-actions/setup-julia/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/setup-julia dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Documenter.yml | 2 +- .github/workflows/Downgrade.yml | 2 +- .github/workflows/Invalidations.yml | 2 +- .github/workflows/benchmark.yml | 2 +- .github/workflows/ci.yml | 2 +- .github/workflows/downstream.yml | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 8f8d674ffa9..1ac862d2c85 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -34,7 +34,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1.9' show-versioninfo: true diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index dd5d8ee7e32..234939a8017 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -64,7 +64,7 @@ jobs: - threaded steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 18048d26be8..b2d34cbc856 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -16,7 +16,7 @@ jobs: if: github.base_ref == github.event.repository.default_branch runs-on: ubuntu-latest steps: - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1' - uses: actions/checkout@v4 diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 4531c3aee0a..15575717f11 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -21,7 +21,7 @@ jobs: - run: | git fetch --tags git branch --create-reflog main origin/main - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a7dfe033a90..27f7af7c528 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -107,7 +107,7 @@ jobs: trixi_test: threaded steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml index b40d5d365c0..83ed6cc79c7 100644 --- a/.github/workflows/downstream.yml +++ b/.github/workflows/downstream.yml @@ -64,7 +64,7 @@ jobs: - TrixiShallowWater.jl steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} From b3bb45745db5803c65ea151a17690b3b4b42dc57 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 3 Apr 2024 08:33:04 +0200 Subject: [PATCH 03/17] 1D Linearized Euler (#1867) * 1D Linearized Euler * remove redundant comment * Fix docu * Update src/equations/linearized_euler_1d.jl * Update examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl * code review * fmt * Apply suggestions from code review Co-authored-by: Andrew Winters * Update src/equations/linearized_euler_1d.jl * Comments * comments * news --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- NEWS.md | 1 + ...r_linearizedeuler_characteristic_system.jl | 112 ++++++++++++++ .../elixir_linearizedeuler_convergence.jl | 59 +++++++ .../elixir_linearizedeuler_gauss_wall.jl | 65 ++++++++ src/Trixi.jl | 2 +- src/equations/equations.jl | 1 + src/equations/linearized_euler_1d.jl | 144 ++++++++++++++++++ test/test_structured_1d.jl | 15 ++ test/test_tree_1d.jl | 3 + test/test_tree_1d_linearizedeuler.jl | 51 +++++++ 10 files changed, 452 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/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 From 876bdeb130099d4d218c817805432f24af1c2285 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 3 Apr 2024 08:50:13 +0200 Subject: [PATCH 04/17] remove duplicate entry from NEWS (#1896) --- NEWS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index e1238562fc9..66dd1e58dc6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,6 @@ for human readability. ## Changes in the v0.7 lifecycle #### Added -- 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. From 01153c503ebc85dd0f5685ed139844af057e7cb5 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 3 Apr 2024 08:51:19 +0200 Subject: [PATCH 05/17] set version to v0.7.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6ff7f29686d..60f340d1306 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.6-pre" +version = "0.7.6" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 90d6d6588cca6a968125c0d1e91a758f566b98a7 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 3 Apr 2024 08:51:32 +0200 Subject: [PATCH 06/17] set development version to v0.7.7-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 60f340d1306..89052859e42 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.6" +version = "0.7.7-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 0f60ebc55189260ec108a07fdbf85d000a7c92d1 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 5 Apr 2024 06:27:08 +0200 Subject: [PATCH 07/17] Math rendering speed of sound docstring linearized Euler Eqs (#1897) * Update linearized_euler_1d.jl * Update linearized_euler_2d.jl --- src/equations/linearized_euler_1d.jl | 2 +- src/equations/linearized_euler_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 22f452219b5..856a53a3d5d 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -24,7 +24,7 @@ Linearized euler equations in one space dimension. The equations are given by 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +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'``. """ diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index d497762bf62..d9d0d44e553 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/src/equations/linearized_euler_2d.jl @@ -29,7 +29,7 @@ Linearized euler equations in two space dimensions. The equations are given by 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +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 LinearizedEulerEquations2D{RealT <: Real} <: From 97aabb7da200c3ed026415e34dfbf1a4afa2623a Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 5 Apr 2024 12:26:18 +0200 Subject: [PATCH 08/17] Fix computation of `max_abs_speeds` for the SWE (#1898) * fix wave speed computation in max_abs_speeds * update cfl values * update test values --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- .../elixir_shallowwater_source_terms.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 2 +- ...xir_shallowwater_quasi_1d_well_balanced.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- .../tree_2d_dgsem/elixir_shallowwater_ec.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- .../elixir_shallowwater_well_balanced_wall.jl | 2 +- .../elixir_shallowwater_ec.jl | 2 +- .../elixir_shallowwater_ec_shockcapturing.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- src/equations/shallow_water_1d.jl | 2 +- src/equations/shallow_water_2d.jl | 2 +- src/equations/shallow_water_quasi_1d.jl | 2 +- test/test_structured_2d.jl | 12 +-- test/test_tree_1d_shallowwater.jl | 18 ++-- test/test_tree_2d_shallowwater.jl | 24 +++--- test/test_unstructured_2d.jl | 82 +++++++++---------- 18 files changed, 84 insertions(+), 80 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl index 48fe37b9996..532fe8dbe7d 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -42,7 +42,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl index a6a56aa807c..09abdf33843 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index a3df37fb966..af0da5d1768 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -77,7 +77,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl index d9f1a52b500..a4f4b0189ba 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl @@ -70,7 +70,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 649e5023f6d..5851530e230 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -73,7 +73,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index bc528ae7756..8221dfebe39 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(dt = 0.2, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl index 13023dfaba2..22043392b2a 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl index f50bd4e4f65..19073b0504a 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl index 9122fb8287d..1f4aa414905 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl index 98408db5a78..c4736e8b9a5 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl @@ -110,7 +110,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl index 6bad3a77f03..6cefca853c1 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -104,7 +104,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index e348ef946b7..7007bea887a 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -551,7 +551,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return (abs(v) + c,) end diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 74a299a51e6..73fae89a0fa 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -927,7 +927,7 @@ end h = waterheight(u, equations) v1, v2 = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return abs(v1) + c, abs(v2) + c end diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 51c360104a7..08620021254 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -248,7 +248,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return (abs(v) + c,) end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 64a1faf05b8..85eb681ec78 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -849,15 +849,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0017285599436729316, - 0.025584610912606776, - 0.028373834961180594, + 0.0017286908591070864, + 0.025585037307655684, + 0.028374244567802766, 6.274146767730866e-5, ], linf=[ - 0.012972309788264802, - 0.108283714215621, - 0.15831585777928936, + 0.012973752001194772, + 0.10829375385832263, + 0.15832858475438094, 0.00018196759554722775, ], tspan=(0.0, 0.05)) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 41ad5c32bbd..8fe3291a938 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -12,10 +12,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[0.244729018751225, 0.8583565222389505, 0.07330427577586297], + l2=[ + 0.24476140682560343, + 0.8587309324660326, + 0.07330427577586297, + ], linf=[ - 2.1635021283528504, - 3.8717508164234453, + 2.1636963952308372, + 3.8737770522883115, 1.7711213427919539, ], tspan=(0.0, 0.25)) @@ -32,13 +36,13 @@ end @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.39464782107209717, - 2.03880864210846, + 0.39472828074570576, + 2.0390687947320076, 4.1623084150546725e-10, ], linf=[ - 0.778905801278281, - 3.2409883402608273, + 0.7793741954662221, + 3.2411927977882096, 7.419800190922032e-10, ], initial_condition=initial_condition_weak_blast_wave, diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 01742644736..9a3ba36c7d4 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -13,15 +13,15 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.991181203601035, - 0.734130029040644, - 0.7447696147162621, + 0.9911802019934329, + 0.7340106828033273, + 0.7446338002084801, 0.5875351036989047, ], linf=[ - 2.0117744577945413, - 2.9962317608172127, - 2.6554999727293653, + 2.0120253138457564, + 2.991158989293406, + 2.6557412817714035, 3.0, ], tspan=(0.0, 0.25)) @@ -280,15 +280,15 @@ end @trixi_testset "elixir_shallowwater_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), l2=[ - 0.13517233723296504, - 0.20010876311162215, - 0.20010876311162223, + 0.1351723240085936, + 0.20010881416550014, + 0.2001088141654999, 2.719538414346464e-7, ], linf=[ - 0.5303607982988336, - 0.5080989745682338, - 0.5080989745682352, + 0.5303608302490757, + 0.5080987791967457, + 0.5080987791967506, 1.1301675764130437e-6, ], tspan=(0.0, 0.25)) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 6814250dd47..35a380e4f24 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -301,15 +301,15 @@ end @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.6106939484178353, - 0.48586236867426724, - 0.48234490854514356, - 0.29467422718511727, + 0.6107326269462766, + 0.48666631722018877, + 0.48309775159067053, + 0.29467422718511704, ], linf=[ - 2.775979948281604, - 3.1721242154451548, - 3.5713448319601393, + 2.776782342826098, + 3.2158378644333707, + 3.652920889487258, 2.052861364219655, ], tspan=(0.0, 0.25)) @@ -408,16 +408,16 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011197623982310795, - 0.04456344888447023, - 0.014317376629669337, - 5.089218476758975e-6, + 0.001118134082248467, + 0.044560486817464634, + 0.01430926600634214, + 5.089218476759981e-6, ], linf=[ - 0.007835284004819698, - 0.3486891284278597, - 0.11242778979399048, - 2.6407324614119432e-5, + 0.007798727223654822, + 0.34782952734839157, + 0.11161614702628064, + 2.6407324614341476e-5, ], tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations @@ -433,15 +433,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.001119678684752799, - 0.015429108794630785, - 0.01708275441241111, - 5.089218476758271e-6, + 0.0011196838135485918, + 0.01542895635133927, + 0.017082803023121197, + 5.089218476759981e-6, ], linf=[ - 0.014299564388827513, - 0.12785126473870534, - 0.17626788561725526, + 0.014299541415654371, + 0.12783948113206955, + 0.17626489583921323, 2.6407324614341476e-5, ], surface_flux=(FluxHydrostaticReconstruction(flux_hll, @@ -461,15 +461,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011196687776346434, - 0.044562672453443995, - 0.014306265289763618, + 0.001118046975499805, + 0.04455969246244461, + 0.014298120235633432, 5.089218476759981e-6, ], linf=[ - 0.007825021762002393, - 0.348550815397918, - 0.1115517935018282, + 0.007776521213640031, + 0.34768318303226353, + 0.11075311228066198, 2.6407324614341476e-5, ], surface_flux=(flux_wintermeyer_etal, @@ -490,15 +490,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011196786847528799, - 0.015429108794631075, - 0.017082754412411742, + 0.0011196838135486059, + 0.015428956351339451, + 0.017082803023120943, 5.089218476759981e-6, ], linf=[ - 0.014299564388830177, - 0.12785126473870667, - 0.17626788561728546, + 0.01429954141565526, + 0.12783948113205668, + 0.176264895839215, 2.6407324614341476e-5, ], surface_flux=(flux_hll, @@ -561,15 +561,15 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec_shockcapturing.jl"), l2=[ - 0.6124656312639043, - 0.504371951785709, - 0.49180896200746366, - 0.29467422718511727, + 0.612551520607341, + 0.5039173660221961, + 0.49136517934903523, + 0.29467422718511704, ], linf=[ - 2.7639232436274392, - 3.3985508653311767, - 3.3330308209196224, + 2.7636771472622197, + 3.236168963021072, + 3.3363936775653826, 2.052861364219655, ], tspan=(0.0, 0.25)) From 6a9cefa16883ea05311c1a5e64bd786769a0b751 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 5 Apr 2024 12:31:49 +0200 Subject: [PATCH 09/17] set version to v0.7.7 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 89052859e42..552f842049d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.7-pre" +version = "0.7.7" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From b4fb1db98d3f3fcc008ce9c6058f59e08ac34a3b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 5 Apr 2024 12:32:05 +0200 Subject: [PATCH 10/17] set development version to v0.7.8-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 552f842049d..5e1a60b7723 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.7" +version = "0.7.8-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From c025873b9bacf26da4f8155f425b3a47bab5f510 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 8 Apr 2024 07:55:49 +0200 Subject: [PATCH 11/17] Boundary integration of viscous forces (#1893) * Viscous force computation * Modularize stress tensor comp for skin friction * typo * try aqua * remove semi from AnalysisSurfaceIntegral * bugfix * two methods for passing aqua tests * rename * test viscous forces * Update src/callbacks_step/analysis_surface_integral_2d.jl Co-authored-by: Andrew Winters * docstring * try to get test to run * Try nicer formatting * revert fmt --------- Co-authored-by: Andrew Winters --- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 18 +- src/Trixi.jl | 3 +- src/callbacks_step/analysis.jl | 30 +- .../analysis_surface_integral_2d.jl | 263 +++++++++++++++--- test/test_parabolic_2d.jl | 29 +- 5 files changed, 296 insertions(+), 47 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 5bfa21c5322..1b485913ab2 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -110,7 +110,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers # Run for a long time to reach a state where forces stabilize up to 3 digits -tspan = (0.0, 1.0) +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) # Callbacks @@ -130,12 +130,26 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, u_inf(equations), l_inf())) +drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + +lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", save_analysis = true, analysis_errors = Symbol[], analysis_integrals = (drag_coefficient, - lift_coefficient)) + lift_coefficient, + drag_coefficient_shear_force, + lift_coefficient_shear_force)) alive_callback = AliveCallback(analysis_interval = analysis_interval) diff --git a/src/Trixi.jl b/src/Trixi.jl index 883f8d66f07..3e7fc5ee519 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, TrivialCallback, AnalysisCallbackCoupled, - AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure + AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, + DragCoefficientShearStress, LiftCoefficientShearStress export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 62048853b53..8f89af755a2 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -9,11 +9,11 @@ # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time """ AnalysisCallback(semi; interval=0, - save_analysis=false, - output_directory="out", - analysis_filename="analysis.dat", - extra_analysis_errors=Symbol[], - extra_analysis_integrals=()) + save_analysis=false, + output_directory="out", + analysis_filename="analysis.dat", + extra_analysis_errors=Symbol[], + extra_analysis_integrals=()) Analyze a numerical solution every `interval` time steps and print the results to the screen. If `save_analysis`, the results are also saved in @@ -634,9 +634,7 @@ pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) # Special analyze for `SemidiscretizationHyperbolicParabolic` such that -# precomputed gradients are available. For now only implemented for the `enstrophy` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +# precomputed gradients are available. function analyze(quantity::typeof(enstrophy), du, u, t, semi::SemidiscretizationHyperbolicParabolic) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @@ -695,3 +693,19 @@ include("analysis_surface_integral_2d.jl") include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") + +# Special analyze for `SemidiscretizationHyperbolicParabolic` such that +# precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces. +# Note that this needs to be included after `analysis_surface_integral_2d.jl` to +# have `VariableViscous` available. +function analyze(quantity::AnalysisSurfaceIntegral{Variable}, + du, u, t, + semi::SemidiscretizationHyperbolicParabolic) where { + Variable <: + VariableViscous} + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + cache_parabolic = semi.cache_parabolic + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, + cache_parabolic) +end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 902fb0b4e85..c5a3d885eb0 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -13,15 +13,19 @@ boundary_symbol_or_boundary_symbols, variable) - This struct is used to compute the surface integral of a quantity of interest `variable` alongside - the boundary/boundaries associated with particular name(s) given in `boundary_symbol` - or `boundary_symbols`. - For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or - drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary - name `:Airfoil` in 2D. +This struct is used to compute the surface integral of a quantity of interest `variable` alongside +the boundary/boundaries associated with particular name(s) given in `boundary_symbol` +or `boundary_symbols`. +For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or +drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary +name `:Airfoil` in 2D. + +- `semi::Semidiscretization`: Passed in to retrieve boundary condition information +- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries + where the quantity of interest is computed +- `variable::Variable`: Quantity of interest, like lift or drag """ -struct AnalysisSurfaceIntegral{Semidiscretization, Variable} - semi::Semidiscretization # passed in to retrieve boundary condition information +struct AnalysisSurfaceIntegral{Variable} indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed variable::Variable # Quantity of interest, like lift or drag @@ -29,8 +33,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} @unpack boundary_symbol_indices = semi.boundary_conditions indices = boundary_symbol_indices[boundary_symbol] - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable) @@ -41,8 +44,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} end sort!(indices) - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end end @@ -50,7 +52,7 @@ struct ForceState{RealT <: Real} psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT - l_inf::RealT + linf::RealT end struct LiftCoefficientPressure{RealT <: Real} @@ -61,13 +63,25 @@ struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end +# Abstract base type used for dispatch of `analyze` for quantities +# requiring gradients of the velocity field. +abstract type VariableViscous end + +struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + +struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + """ - LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + LiftCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the lift coefficient ```math C_{L,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_L \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -76,25 +90,25 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function LiftCoefficientPressure(aoa, rhoinf, uinf, linf) # psi_lift is the normal unit vector to the freestream direction. # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) # leads to positive lift coefficients for positive angles of attack for airfoils. # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same # value, but with the opposite sign. psi_lift = (-sin(aoa), cos(aoa)) - return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, l_inf)) + return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf)) end """ - DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + DragCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the drag coefficient ```math C_{D,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_D \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -103,28 +117,140 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, linf) + # `psi_drag` is the unit vector tangent to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf)) +end + +""" + LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the lift coefficient +```math +C_{L,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_L \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + # psi_lift is the normal unit vector to the freestream direction. + # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) + # leads to negative lift coefficients for airfoils. + # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same + # value, but with the opposite sign. + psi_lift = (-sin(aoa), cos(aoa)) + return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf)) +end + +""" + DragCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the drag coefficient +```math +C_{D,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_D \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function DragCoefficientShearStress(aoa, rhoinf, uinf, linf) # `psi_drag` is the unit vector tangent to the freestream direction psi_drag = (cos(aoa), sin(aoa)) - return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf)) + return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf)) end function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = lift_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) end function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = drag_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) +end + +# Compute the three components of the symmetric viscous stress tensor +# (tau_11, tau_12, tau_22) based on the gradients of the velocity field. +# This is required for drag and lift coefficients based on shear stress, +# as well as for the non-integrated quantities such as +# skin friction coefficient (to be added). +function viscous_stress_tensor(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients_1, + equations_parabolic) + _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients_2, + equations_parabolic) + + # Components of viscous stress tensor + # (4/3 * (v1)_x - 2/3 * (v2)_y) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy + # ((v1)_y + (v2)_x) + # stress tensor is symmetric + tau_12 = dv1dy + dv2dx # = tau_21 + # (4/3 * (v2)_y - 2/3 * (v1)_x) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx + + mu = dynamic_viscosity(u, equations_parabolic) + + return mu .* (tau_11, tau_12, tau_22) +end + +function viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + # Normalize normal direction, should point *into* the fluid => *(-1) + n_normal = -normal_direction / norm(normal_direction) + + tau_11, tau_12, tau_22 = viscous_stress_tensor(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + + # Viscous stress vector: Stress tensor * normal vector + visc_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2] + visc_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2] + + return (visc_stress_vector_1, visc_stress_vector_2) +end + +function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state + return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / + (0.5 * rhoinf * uinf^2 * linf) +end + +function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state + return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / + (0.5 * rhoinf * uinf^2 * linf) end function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @@ -169,21 +295,88 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, + du, u, t, mesh::P4estMesh{2}, + equations, equations_parabolic, + dg::DGSEM, cache, + cache_parabolic) where {Variable <: VariableViscous} + @unpack boundaries = cache + @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements + @unpack weights = dg.basis + + @unpack indices, variable = surface_variable + + # Additions for parabolic + @unpack viscous_container = cache_parabolic + @unpack gradients = viscous_container + + gradients_x, gradients_y = gradients + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + for boundary in indices + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in index_range + u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, + boundary) + # Extract normal direction at nodes which points from the elements outwards, + # i.e., *into* the structure. + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, + element) + + # L2 norm of normal direction (contravariant_vector) is the surface element + dS = weights[node_index] * norm(normal_direction) + + gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node, + j_node, element) + gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node, + j_node, element) + + # Integral over whole boundary surface + surface_integral += variable(u_node, normal_direction, equations_parabolic, + gradients_1, gradients_2) * dS + + i_node += i_node_step + j_node += j_node_step + end + end + return surface_integral +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end end # muladd diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index edbdfe5a84f..a18dfb353b7 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -719,7 +719,10 @@ end 1.199362305026636, 0.9077214424040279, 5.666071182328691], tspan=(0.0, 0.001), - initial_refinement_level=0) + initial_refinement_level=0, + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 10_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -728,6 +731,30 @@ end du_ode = similar(u_ode) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + + @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13) + @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13) + + @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13) + @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13) end end From c6fc9c591726e6d6d51605f46d879408202d625d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 9 Apr 2024 13:42:44 +0200 Subject: [PATCH 12/17] Use `eachindex` instead of `1:length`, compare against `nothing` with `===` (#1899) * Use `eachindex` instead of `1:length` and compare against nothing with `===` * revert some places * revert for solvers * revert * revert * revert * fmt * fmt * revert --- .../src/files/non_periodic_boundaries.jl | 2 +- .../src/files/scalar_linear_advection_1d.jl | 4 ++-- .../elixir_advection_amr_solution_independent.jl | 2 +- .../elixir_advection_amr_solution_independent.jl | 2 +- .../elixir_advection_amr_solution_independent.jl | 2 +- src/callbacks_step/amr.jl | 8 ++++---- .../euler_acoustics_coupling_dg2d.jl | 2 +- src/callbacks_step/time_series_dg_tree.jl | 16 ++++++++-------- .../time_series_dg_unstructured.jl | 10 +++++----- .../semidiscretization_euler_acoustics.jl | 4 ++-- src/visualization/utilities.jl | 4 ++-- test/test_unit.jl | 2 +- utils/trixi2txt.jl | 8 ++++---- 13 files changed, 33 insertions(+), 33 deletions(-) diff --git a/docs/literate/src/files/non_periodic_boundaries.jl b/docs/literate/src/files/non_periodic_boundaries.jl index 3b238ad4533..8f0e320dfdc 100644 --- a/docs/literate/src/files/non_periodic_boundaries.jl +++ b/docs/literate/src/files/non_periodic_boundaries.jl @@ -99,7 +99,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), save_everystep=false, saveat=visnodes, callback=callbacks); using Plots -@gif for step in 1:length(sol.u) +@gif for step in eachindex(sol.u) plot(sol.u[step], semi, ylim=(-1.5, 1.5), legend=true, label="approximation", title="time t=$(round(sol.t[step], digits=5))") scatter!([0.0], [sum(boundary_condition(SVector(0.0), sol.t[step], equations))], label="boundary condition") end diff --git a/docs/literate/src/files/scalar_linear_advection_1d.jl b/docs/literate/src/files/scalar_linear_advection_1d.jl index 9b48f29d341..c7d55e26d2a 100644 --- a/docs/literate/src/files/scalar_linear_advection_1d.jl +++ b/docs/literate/src/files/scalar_linear_advection_1d.jl @@ -120,7 +120,7 @@ integral = sum(nodes.^3 .* weights) x = Matrix{Float64}(undef, length(nodes), n_elements) for element in 1:n_elements x_l = coordinates_min + (element - 1) * dx + dx/2 - for i in 1:length(nodes) + for i in eachindex(nodes) ξ = nodes[i] # nodes in [-1, 1] x[i, element] = x_l + dx/2 * ξ end @@ -417,7 +417,7 @@ dx = (coordinates_max - coordinates_min) / n_elements # length of one element x = Matrix{Float64}(undef, length(nodes), n_elements) for element in 1:n_elements x_l = -1 + (element - 1) * dx + dx/2 - for i in 1:length(nodes) # basis points in [-1, 1] + for i in eachindex(nodes) # basis points in [-1, 1] ξ = nodes[i] x[i, element] = x_l + dx/2 * ξ end diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl index 5a2537be4e6..a1ddc6a314f 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -33,7 +33,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 #Iterate over all elements - for element in 1:length(alpha) + for element in eachindex(alpha) # Calculate periodic distance between cell and center. # This requires an uncurved mesh! coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] + diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl index 0589e76a6a9..1ed08e1961b 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -32,7 +32,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 # Iterate over all elements. - for element in 1:length(alpha) + for element in eachindex(alpha) # Calculate periodic distance between cell and center. # This requires an uncurved mesh! coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] + diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl index 03a213689ec..c7412660b0c 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -31,7 +31,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 #Iterate over all elements - for element in 1:length(alpha) + for element in eachindex(alpha) #Calculate periodic distance between cell and center. cell_id = cache.elements.cell_ids[element] coordinates = mesh.tree.coordinates[1:2, cell_id] diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 1ab65a3553e..45f03fba8fe 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -243,7 +243,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @unpack to_refine, to_coarsen = amr_callback.amr_cache empty!(to_refine) empty!(to_coarsen) - for element in 1:length(lambda) + for element in eachindex(lambda) controller_value = lambda[element] if controller_value > 0 push!(to_refine, leaf_cell_ids[element]) @@ -307,7 +307,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, end # Extract only those parent cells for which all children should be coarsened - to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] + to_coarsen = collect(eachindex(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] # Finally, coarsen mesh coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh.tree, @@ -395,7 +395,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @unpack to_refine, to_coarsen = amr_callback.amr_cache empty!(to_refine) empty!(to_coarsen) - for element in 1:length(lambda) + for element in eachindex(lambda) controller_value = lambda[element] if controller_value > 0 push!(to_refine, leaf_cell_ids[element]) @@ -456,7 +456,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, end # Extract only those parent cells for which all children should be coarsened - to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] + to_coarsen = collect(eachindex(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] # Finally, coarsen mesh coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh.tree, diff --git a/src/callbacks_step/euler_acoustics_coupling_dg2d.jl b/src/callbacks_step/euler_acoustics_coupling_dg2d.jl index 16fac4f2d8d..8a8bb893dcd 100644 --- a/src/callbacks_step/euler_acoustics_coupling_dg2d.jl +++ b/src/callbacks_step/euler_acoustics_coupling_dg2d.jl @@ -12,7 +12,7 @@ function calc_acoustic_sources!(acoustic_source_terms, u_euler, u_acoustics, dg::DGSEM, cache) acoustic_source_terms .= zero(eltype(acoustic_source_terms)) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) diff --git a/src/callbacks_step/time_series_dg_tree.jl b/src/callbacks_step/time_series_dg_tree.jl index 37d4e6ea705..0af1688a8ed 100644 --- a/src/callbacks_step/time_series_dg_tree.jl +++ b/src/callbacks_step/time_series_dg_tree.jl @@ -25,7 +25,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, cell_id = cell_ids[element] # Iterate over coordinates - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Skip coordinates for which an element has already been found if element_ids[index] > 0 continue @@ -63,7 +63,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, wbary = barycentric_weights(nodes) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Construct point x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) @@ -94,7 +94,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -108,7 +108,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index]) end @@ -126,7 +126,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -140,7 +140,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index]) @@ -159,7 +159,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -173,7 +173,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, k, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index] diff --git a/src/callbacks_step/time_series_dg_unstructured.jl b/src/callbacks_step/time_series_dg_unstructured.jl index f6d1bb48f24..85427f1273a 100644 --- a/src/callbacks_step/time_series_dg_unstructured.jl +++ b/src/callbacks_step/time_series_dg_unstructured.jl @@ -31,7 +31,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Iterate over coordinates distances = zeros(eltype(mesh.corners), mesh.n_elements) indices = zeros(Int, mesh.n_elements, 2) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Grab the current point for which the element needs found point = SVector(coordinates[1, index], coordinates[2, index]) @@ -77,7 +77,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Loop through all the element candidates until we find a vector from the barycenter # to the surface that points in the same direction as the current `point` vector. # This then gives us the correct element. - for element in 1:length(candidates) + for element in eachindex(candidates) bary_center = SVector(bary_centers[1, candidates[element]], bary_centers[2, candidates[element]]) # Vector pointing from the barycenter toward the minimal `surface_point` @@ -153,7 +153,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, # Helper array for a straight-sided quadrilateral element corners = zeros(eltype(mesh.corners), 4, 2) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Construct point x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) @@ -280,7 +280,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -294,7 +294,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index]) diff --git a/src/semidiscretization/semidiscretization_euler_acoustics.jl b/src/semidiscretization/semidiscretization_euler_acoustics.jl index 173523ff892..286315fb960 100644 --- a/src/semidiscretization/semidiscretization_euler_acoustics.jl +++ b/src/semidiscretization/semidiscretization_euler_acoustics.jl @@ -103,7 +103,7 @@ function precompute_weights(source_region, weights, coupled_element_ids, equatio (nnodes(dg), nnodes(dg), length(coupled_element_ids))) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) x = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, @@ -197,7 +197,7 @@ function add_acoustic_source_terms!(du_acoustics, acoustic_source_terms, source_ coupled_element_ids, mesh::TreeMesh{2}, equations, dg::DGSEM, cache) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index 05457395ac0..c1108128c92 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -495,7 +495,7 @@ function cell2node(cell_centered_data) resolution_in, _ = size(first(cell_centered_data)) resolution_out = resolution_in + 1 node_centered_data = [Matrix{Float64}(undef, resolution_out, resolution_out) - for _ in 1:length(cell_centered_data)] + for _ in eachindex(cell_centered_data)] for (cell_data, node_data) in zip(cell_centered_data, node_centered_data) # Fill center with original data @@ -1545,7 +1545,7 @@ end # Create an axis. function axis_curve(nodes_x, nodes_y, nodes_z, slice, point, n_points) - if n_points == nothing + if n_points === nothing n_points = 64 end dimensions = length(point) diff --git a/test/test_unit.jl b/test/test_unit.jl index 79950f58d59..9a4ae2ede13 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -284,7 +284,7 @@ end function MyContainer(data, capacity) c = MyContainer(Vector{Int}(undef, capacity + 1), capacity, length(data), capacity + 1) - c.data[1:length(data)] .= data + c.data[eachindex(data)] .= data return c end MyContainer(data::AbstractArray) = MyContainer(data, length(data)) diff --git a/utils/trixi2txt.jl b/utils/trixi2txt.jl index 12a3d46760e..52ee904d2f6 100644 --- a/utils/trixi2txt.jl +++ b/utils/trixi2txt.jl @@ -86,7 +86,7 @@ function trixi2txt(filename::AbstractString...; "maximum supported level $max_supported_level") end max_available_nodes_per_finest_element = 2^(max_supported_level - max_level) - if nvisnodes == nothing + if nvisnodes === nothing max_nvisnodes = 2 * n_nodes elseif nvisnodes == 0 max_nvisnodes = n_nodes @@ -137,9 +137,9 @@ function trixi2txt(filename::AbstractString...; println(io) # Data - for idx in 1:length(xs) + for idx in eachindex(xs) @printf(io, "%+10.8e", xs[idx]) - for variable_id in 1:length(variables) + for variable_id in eachindex(variables) @printf(io, " %+10.8e ", node_centered_data[idx, variable_id]) end println(io) @@ -199,7 +199,7 @@ function read_meshfile(filename::String) # Extract leaf cells (= cells to be plotted) and contract all other arrays accordingly leaf_cells = similar(levels) n_cells = 0 - for cell_id in 1:length(levels) + for cell_id in eachindex(levels) if sum(child_ids[:, cell_id]) > 0 continue end From efbd5a7f740d7b81fc2b2a8c1d4662988c0d8f83 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Tue, 9 Apr 2024 16:51:54 +0200 Subject: [PATCH 13/17] Add preprocessing step to automate numbering in tutorial introduction (#1901) * Add preprocessing of tutorial introduction file * Rename variables * implement suggestions --- docs/literate/make.jl | 12 ++++++- docs/literate/src/files/index.jl | 60 +++++++++++++++++++++----------- 2 files changed, 51 insertions(+), 21 deletions(-) diff --git a/docs/literate/make.jl b/docs/literate/make.jl index 84e4fbdced6..262a236971c 100644 --- a/docs/literate/make.jl +++ b/docs/literate/make.jl @@ -75,7 +75,17 @@ function create_tutorials(files) end # Generate markdown file for introduction page - Literate.markdown(joinpath(repo_src, "index.jl"), pages_dir; name="introduction") + # Preprocessing introduction file: Generate consecutive tutorial numbers by replacing + # each occurrence of `{index}` with an integer incremented by 1, starting at 1. + function preprocess_introduction(content) + counter = 1 + while occursin("{index}", content) + content = replace(content, "{index}" => "$counter", count = 1) + counter += 1 + end + return content + end + Literate.markdown(joinpath(repo_src, "index.jl"), pages_dir; name="introduction", preprocess=preprocess_introduction) # Navigation system for makedocs pages = Any["Introduction" => "tutorials/introduction.md",] diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index 6d5800c3b66..5605803db22 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -14,14 +14,16 @@ # There are tutorials for the following topics: -# ### [1 First steps in Trixi.jl](@ref getting_started) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} First steps in Trixi.jl](@ref getting_started) #- # This tutorial provides guidance for getting started with Trixi.jl, and Julia as well. It outlines # the installation procedures for both Julia and Trixi.jl, the execution of Trixi.jl elixirs, the # fundamental structure of a Trixi.jl setup, the visualization of results, and the development # process for Trixi.jl. -# ### [2 Behind the scenes of a simulation setup](@ref behind_the_scenes_simulation_setup) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Behind the scenes of a simulation setup](@ref behind_the_scenes_simulation_setup) #- # This tutorial will guide you through a simple Trixi.jl setup ("elixir"), giving an overview of # what happens in the background during the initialization of a simulation. While the setup @@ -30,20 +32,23 @@ # the more fundamental, *technical* concepts that are applicable to a variety of # (also more complex) configurations.s -# ### [3 Introduction to DG methods](@ref scalar_linear_advection_1d) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Introduction to DG methods](@ref scalar_linear_advection_1d) #- # This tutorial gives an introduction to discontinuous Galerkin (DG) methods with the example of the # scalar linear advection equation in 1D. Starting with some theoretical explanations, we first implement # a raw version of a discontinuous Galerkin spectral element method (DGSEM). Then, we will show how # to use features of Trixi.jl to achieve the same result. -# ### [4 DGSEM with flux differencing](@ref DGSEM_FluxDiff) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} DGSEM with flux differencing](@ref DGSEM_FluxDiff) #- # To improve stability often the flux differencing formulation of the DGSEM (split form) is used. # We want to present the idea and formulation on a basic 1D level. Then, we show how this formulation # can be implemented in Trixi.jl and analyse entropy conservation for two different flux combinations. -# ### [5 Shock capturing with flux differencing and stage limiter](@ref shock_capturing) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Shock capturing with flux differencing and stage limiter](@ref shock_capturing) #- # Using the flux differencing formulation, a simple procedure to capture shocks is a hybrid blending # of a high-order DG method and a low-order subcell finite volume (FV) method. We present the idea on a @@ -51,20 +56,23 @@ # explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler # equations. -# ### [6 Non-periodic boundary conditions](@ref non_periodic_boundaries) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Non-periodic boundary conditions](@ref non_periodic_boundaries) #- # Thus far, all examples used periodic boundaries. In Trixi.jl, you can also set up a simulation with # non-periodic boundaries. This tutorial presents the implementation of the classical Dirichlet # boundary condition with a following example. Then, other non-periodic boundaries are mentioned. -# ### [7 DG schemes via `DGMulti` solver](@ref DGMulti_1) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} DG schemes via `DGMulti` solver](@ref DGMulti_1) #- # This tutorial is about the more general DG solver [`DGMulti`](@ref), introduced [here](@ref DGMulti). # We are showing some examples for this solver, for instance with discretization nodes by Gauss or # triangular elements. Moreover, we present a simple way to include pre-defined triangulate meshes for # non-Cartesian domains using the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl). -# ### [8 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2) #- # Supplementary to the previous tutorial about DG schemes via the `DGMulti` solver we now present # the possibility for `DGMulti` to use other SBP schemes via the package @@ -72,7 +80,8 @@ # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin # (CGSEM) method. -# ### [9 Upwind FD SBP schemes](@ref upwind_fdsbp) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Upwind FD SBP schemes](@ref upwind_fdsbp) #- # General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but # also with a general `DG` solver. In particular, upwind finite difference SBP @@ -80,42 +89,49 @@ # schemes in the `DGMulti` framework, the interface is based on the package # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl). -# ### [10 Adding a new scalar conservation law](@ref adding_new_scalar_equations) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Adding a new scalar conservation law](@ref adding_new_scalar_equations) #- # This tutorial explains how to add a new physics model using the example of the cubic conservation # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then, # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented # and the ODE problem is solved by OrdinaryDiffEq's `solve` method. -# ### [11 Adding a non-conservative equation](@ref adding_nonconservative_equation) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Adding a non-conservative equation](@ref adding_nonconservative_equation) #- # In this part, another physics model is implemented, the nonconservative linear advection equation. # We run two different simulations with different levels of refinement and compare the resulting errors. -# ### [12 Parabolic terms](@ref parabolic_terms) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Parabolic terms](@ref parabolic_terms) #- # This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g., # to solve the advection-diffusion equation. -# ### [13 Adding new parabolic terms](@ref adding_new_parabolic_terms) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Adding new parabolic terms](@ref adding_new_parabolic_terms) #- # This tutorial describes how new parabolic terms can be implemented using Trixi.jl. -# ### [14 Adaptive mesh refinement](@ref adaptive_mesh_refinement) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Adaptive mesh refinement](@ref adaptive_mesh_refinement) #- # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while # not wasting resources for less interesting parts of the domain. This leads to much more efficient # simulations. This tutorial presents the implementation strategy of AMR in Trixi.jl, including the use of # different indicators and controllers. -# ### [15 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Structured mesh with curvilinear mapping](@ref structured_mesh_mapping) #- # In this tutorial, the use of Trixi.jl's structured curved mesh type [`StructuredMesh`](@ref) is explained. # We present the two basic option to initialize such a mesh. First, the curved domain boundaries # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is # defined by passing the transformation mapping. -# ### [16 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial) #- # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref) # functionality of Trixi.jl. This begins by running and visualizing an available unstructured @@ -124,26 +140,30 @@ # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh. # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`. -# ### [17 P4est mesh from gmsh](@ref p4est_from_gmsh) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} P4est mesh from gmsh](@ref p4est_from_gmsh) #- # This tutorial describes how to obtain a [`P4estMesh`](@ref) from an existing mesh generated # by [`gmsh`](https://gmsh.info/) or any other meshing software that can export to the Abaqus # input `.inp` format. The tutorial demonstrates how edges/faces can be associated with boundary conditions based on the physical nodesets. -# ### [18 Explicit time stepping](@ref time_stepping) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Explicit time stepping](@ref time_stepping) #- # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl). # It explains how to use their algorithms and presents two types of time step choices - with error-based # and CFL-based adaptive step size control. -# ### [19 Differentiable programming](@ref differentiable_programming) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Differentiable programming](@ref differentiable_programming) #- # This part deals with some basic differentiable programming topics. For example, a Jacobian, its # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for # a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl # at the end. -# ### [20 Custom semidiscretization](@ref custom_semidiscretization) +#src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing +# ### [{index} Custom semidiscretization](@ref custom_semidiscretization) #- # This tutorial describes the [semidiscretiations](@ref overview-semidiscretizations) of Trixi.jl # and explains how to extend them for custom tasks. From b1a84a63298966b67491dacac78e785e0dbfb36c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 10 Apr 2024 09:45:04 +0200 Subject: [PATCH 14/17] add x coordinates and time as input arguments to the variable function to be integrated along a boundary (#1902) --- .../analysis_surface_integral_2d.jl | 42 +++++++++++++------ 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index c5a3d885eb0..7ae259e5285 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -63,7 +63,7 @@ struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end -# Abstract base type used for dispatch of `analyze` for quantities +# Abstract base type used for dispatch of `analyze` for quantities # requiring gradients of the velocity field. abstract type VariableViscous end @@ -175,7 +175,8 @@ function DragCoefficientShearStress(aoa, rhoinf, uinf, linf) return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf)) end -function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) +function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, x, t, + equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector @@ -183,7 +184,8 @@ function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equati return p * n / (0.5 * rhoinf * uinf^2 * linf) end -function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) +function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, x, t, + equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector @@ -193,8 +195,8 @@ end # Compute the three components of the symmetric viscous stress tensor # (tau_11, tau_12, tau_22) based on the gradients of the velocity field. -# This is required for drag and lift coefficients based on shear stress, -# as well as for the non-integrated quantities such as +# This is required for drag and lift coefficients based on shear stress, +# as well as for the non-integrated quantities such as # skin friction coefficient (to be added). function viscous_stress_tensor(u, normal_direction, equations_parabolic, gradients_1, gradients_2) @@ -233,7 +235,7 @@ function viscous_stress_vector(u, normal_direction, equations_parabolic, return (visc_stress_vector_1, visc_stress_vector_2) end -function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, +function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x, t, equations_parabolic, gradients_1, gradients_2) visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, @@ -243,7 +245,7 @@ function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, (0.5 * rhoinf * uinf^2 * linf) end -function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, +function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, t, equations_parabolic, gradients_1, gradients_2) visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, @@ -283,10 +285,18 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, i_node, j_node, element) + # Coordinates at a boundary node + x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, + element) + # L2 norm of normal direction (contravariant_vector) is the surface element dS = weights[node_index] * norm(normal_direction) - # Integral over entire boundary surface - surface_integral += variable(u_node, normal_direction, equations) * dS + + # Integral over entire boundary surface. Note, it is assumed that the + # `normal_direction` is normalized to be a normal vector within the + # function `variable` and the division of the normal scaling factor + # `norm(normal_direction)` is then accounted for with the `dS` quantity. + surface_integral += variable(u_node, normal_direction, x, t, equations) * dS i_node += i_node_step j_node += j_node_step @@ -328,11 +338,15 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, boundary) # Extract normal direction at nodes which points from the elements outwards, - # i.e., *into* the structure. + # i.e., *into* the structure. normal_direction = get_normal_direction(direction, contravariant_vectors, i_node, j_node, element) + # Coordinates at a boundary node + x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, + element) + # L2 norm of normal direction (contravariant_vector) is the surface element dS = weights[node_index] * norm(normal_direction) @@ -341,8 +355,12 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node, j_node, element) - # Integral over whole boundary surface - surface_integral += variable(u_node, normal_direction, equations_parabolic, + # Integral over whole boundary surface. Note, it is assumed that the + # `normal_direction` is normalized to be a normal vector within the + # function `variable` and the division of the normal scaling factor + # `norm(normal_direction)` is then accounted for with the `dS` quantity. + surface_integral += variable(u_node, normal_direction, x, t, + equations_parabolic, gradients_1, gradients_2) * dS i_node += i_node_step From 3cd1266d8ef2d470c1aa3ec802696727e9293ec7 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 16 Apr 2024 06:21:41 +0200 Subject: [PATCH 15/17] LinearizedEuler3D (#1903) * LinearizedEuler3D * fmt * Apply suggestions from code review Co-authored-by: Lars Christmann * nicer IC * test vals * fmt * unit tests for full coverage --------- Co-authored-by: Lars Christmann --- .../elixir_linearizedeuler_convergence.jl | 64 +++++ .../elixir_linearizedeuler_gauss_wall.jl | 65 +++++ src/Trixi.jl | 2 +- src/equations/equations.jl | 1 + src/equations/linearized_euler_1d.jl | 2 +- src/equations/linearized_euler_2d.jl | 2 +- src/equations/linearized_euler_3d.jl | 254 ++++++++++++++++++ test/test_p4est_3d.jl | 22 ++ test/test_tree_3d_linearizedeuler.jl | 35 +++ test/test_tree_3d_part2.jl | 3 + test/test_unit.jl | 54 ++++ 11 files changed, 501 insertions(+), 3 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl create mode 100644 examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl create mode 100644 src/equations/linearized_euler_3d.jl create mode 100644 test/test_tree_3d_linearizedeuler.jl diff --git a/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl b/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl new file mode 100644 index 00000000000..fd066ef8955 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_linearizedeuler_convergence.jl @@ -0,0 +1,64 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations3D(v_mean_global = (0.0, 0.0, 0.0), c_mean_global = 1.0, + rho_mean_global = 1.0) + +initial_condition = initial_condition_convergence_test + +solver = DGSEM(polydeg = 3, surface_flux = flux_hll) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# `initial_refinement_level` is provided here to allow for a +# convenient convergence test, see +# https://trixi-framework.github.io/Trixi.jl/stable/#Performing-a-convergence-analysis +trees_per_dimension = (4, 4, 4) +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 0.2 +tspan = (0.0, 0.2) +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_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl new file mode 100644 index 00000000000..2eb8ef822a3 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -0,0 +1,65 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations3D(v_mean_global = (0.5, 0.5, 0.5), c_mean_global = 1.0, + rho_mean_global = 1.0) + +solver = DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (90.0, 90.0, 90.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000, + periodicity = false) + +# Initialize density and pressure perturbation with a Gaussian bump +# that splits into radial waves which are advected with v - c and v + c. +function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations3D) + v1_prime = 0.0 + v2_prime = 0.0 + v3_prime = 0.0 + rho_prime = p_prime = 2 * exp(-((x[1] - 45)^2 + (x[2] - 45)^2) / 25) + return SVector(rho_prime, v1_prime, v2_prime, v3_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. + +# At t = 30, the wave moving with v + c crashes into the wall +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.9) + +# 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 3e7fc5ee519..5511f7e02e2 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -160,7 +160,7 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations1D, LinearizedEulerEquations2D, + LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D, PolytropicEulerEquations2D, TrafficFlowLWREquations1D diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ef4dac8b1a5..69a54f7faf8 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -503,6 +503,7 @@ abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("linearized_euler_1d.jl") include("linearized_euler_2d.jl") +include("linearized_euler_3d.jl") abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <: AbstractEquations{NDIMS, NVARS} end diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 856a53a3d5d..19a3bdcb3bd 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -8,7 +8,7 @@ @doc raw""" LinearizedEulerEquations1D(v_mean_global, c_mean_global, rho_mean_global) -Linearized euler equations in one space dimension. The equations are given by +Linearized Euler equations in one space dimension. The equations are given by ```math \partial_t \begin{pmatrix} diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index d9d0d44e553..3df3093069d 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/src/equations/linearized_euler_2d.jl @@ -8,7 +8,7 @@ @doc raw""" LinearizedEulerEquations2D(v_mean_global, c_mean_global, rho_mean_global) -Linearized euler equations in two space dimensions. The equations are given by +Linearized Euler equations in two space dimensions. The equations are given by ```math \partial_t \begin{pmatrix} diff --git a/src/equations/linearized_euler_3d.jl b/src/equations/linearized_euler_3d.jl new file mode 100644 index 00000000000..ab5f9863db8 --- /dev/null +++ b/src/equations/linearized_euler_3d.jl @@ -0,0 +1,254 @@ +# 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""" + LinearizedEulerEquations3D(v_mean_global, c_mean_global, rho_mean_global) + +Linearized Euler equations in three space dimensions. The equations are given by +```math +\partial_t +\begin{pmatrix} + \rho' \\ v_1' \\ v_2' \\ v_3' \\ 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} v_2' \\ + \bar{v_1} v_3' \\ + \bar{v_1} p' + c^2 \bar{\rho} v_1' +\end{pmatrix} ++ +\partial_y +\begin{pmatrix} + \bar{\rho} v_2' + \bar{v_2} \rho ' \\ + \bar{v_2} v_1' \\ + \bar{v_2} v_2' + \frac{p'}{\bar{\rho}} \\ + \bar{v_2} v_3' \\ + \bar{v_2} p' + c^2 \bar{\rho} v_2' +\end{pmatrix} ++ +\partial_z +\begin{pmatrix} + \bar{\rho} v_3' + \bar{v_3} \rho ' \\ + \bar{v_3} v_1' \\ + \bar{v_3} v_2' \\ + \bar{v_3} v_3' + \frac{p'}{\bar{\rho}} \\ + \bar{v_3} p' + c^2 \bar{\rho} v_3' +\end{pmatrix} += +\begin{pmatrix} + 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, v_3')``, the pressure ``p'`` and the density ``\rho'``. +""" +struct LinearizedEulerEquations3D{RealT <: Real} <: + AbstractLinearizedEulerEquations{3, 5} + v_mean_global::SVector{3, RealT} + c_mean_global::RealT + rho_mean_global::RealT +end + +function LinearizedEulerEquations3D(v_mean_global::NTuple{3, <: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 LinearizedEulerEquations3D(SVector(v_mean_global), c_mean_global, + rho_mean_global) +end + +function LinearizedEulerEquations3D(; v_mean_global::NTuple{3, <:Real}, + c_mean_global::Real, rho_mean_global::Real) + return LinearizedEulerEquations3D(v_mean_global, c_mean_global, + rho_mean_global) +end + +function varnames(::typeof(cons2cons), ::LinearizedEulerEquations3D) + ("rho_prime", "v1_prime", "v2_prime", "v3_prime", "p_prime") +end +function varnames(::typeof(cons2prim), ::LinearizedEulerEquations3D) + ("rho_prime", "v1_prime", "v2_prime", "v3_prime", "p_prime") +end + +""" + initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations3D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations3D) + rho_prime = -cospi(2 * t) * (sinpi(2 * x[1]) + sinpi(2 * x[2]) + sinpi(2 * x[3])) + v1_prime = sinpi(2 * t) * cospi(2 * x[1]) + v2_prime = sinpi(2 * t) * cospi(2 * x[2]) + v3_prime = sinpi(2 * t) * cospi(2 * x[3]) + p_prime = rho_prime + + return SVector(rho_prime, v1_prime, v2_prime, v3_prime, p_prime) +end + +""" + boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function, + equations::LinearizedEulerEquations3D) + +Boundary conditions for a solid wall. +""" +function boundary_condition_wall(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::LinearizedEulerEquations3D) + # 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. + # Similarly, for boundaries in the -y/+y or -z/+z direction, we multiply the + # velocity in the y or z direction by -1 + if direction in (1, 2) # x direction + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4], + u_inner[5]) + elseif direction in (3, 4) # y direction + u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4], + u_inner[5]) + else # z direction = (5, 6) + u_boundary = SVector(u_inner[1], u_inner[2], u_inner[3], -u_inner[4], + u_inner[5]) + end + + # Calculate boundary flux depending on the orientation of the boundary + # Odd directions are in negative coordinate direction, + # even directions are in positive coordinate direction. + 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 3D flux for a single point +@inline function flux(u, orientation::Integer, equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, v2_prime, v3_prime, p_prime = u + if orientation == 1 + f1 = v_mean_global[1] * rho_prime + rho_mean_global * v1_prime + f2 = v_mean_global[1] * v1_prime + p_prime / rho_mean_global + f3 = v_mean_global[1] * v2_prime + f4 = v_mean_global[1] * v3_prime + f5 = v_mean_global[1] * p_prime + c_mean_global^2 * rho_mean_global * v1_prime + elseif orientation == 2 + f1 = v_mean_global[2] * rho_prime + rho_mean_global * v2_prime + f2 = v_mean_global[2] * v1_prime + f3 = v_mean_global[2] * v2_prime + p_prime / rho_mean_global + f4 = v_mean_global[2] * v3_prime + f5 = v_mean_global[2] * p_prime + c_mean_global^2 * rho_mean_global * v2_prime + else # orientation == 3 + f1 = v_mean_global[3] * rho_prime + rho_mean_global * v3_prime + f2 = v_mean_global[3] * v1_prime + f3 = v_mean_global[3] * v2_prime + f4 = v_mean_global[3] * v3_prime + p_prime / rho_mean_global + f5 = v_mean_global[3] * p_prime + c_mean_global^2 * rho_mean_global * v3_prime + end + + return SVector(f1, f2, f3, f4, f5) +end + +# Calculate 3D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, v2_prime, v3_prime, p_prime = u + + v_mean_normal = v_mean_global[1] * normal_direction[1] + + v_mean_global[2] * normal_direction[2] + + v_mean_global[3] * normal_direction[3] + v_prime_normal = v1_prime * normal_direction[1] + v2_prime * normal_direction[2] + + v3_prime * normal_direction[3] + + f1 = v_mean_normal * rho_prime + rho_mean_global * v_prime_normal + f2 = v_mean_normal * v1_prime + normal_direction[1] * p_prime / rho_mean_global + f3 = v_mean_normal * v2_prime + normal_direction[2] * p_prime / rho_mean_global + f4 = v_mean_normal * v3_prime + normal_direction[3] * p_prime / rho_mean_global + f5 = v_mean_normal * p_prime + c_mean_global^2 * rho_mean_global * v_prime_normal + + return SVector(f1, f2, f3, f4, f5) +end + +@inline have_constant_speed(::LinearizedEulerEquations3D) = True() + +@inline function max_abs_speeds(equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global[1]) + c_mean_global, abs(v_mean_global[2]) + c_mean_global, + abs(v_mean_global[3]) + c_mean_global +end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + if orientation == 1 + return abs(v_mean_global[1]) + c_mean_global + elseif orientation == 2 + return abs(v_mean_global[2]) + c_mean_global + else # orientation == 3 + return abs(v_mean_global[3]) + c_mean_global + end +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + v_mean_normal = normal_direction[1] * v_mean_global[1] + + normal_direction[2] * v_mean_global[2] + + normal_direction[3] * v_mean_global[3] + return abs(v_mean_normal) + c_mean_global * norm(normal_direction) +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::LinearizedEulerEquations3D) + min_max_speed_davis(u_ll, u_rr, orientation, equations) +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + min_max_speed_davis(u_ll, u_rr, normal_direction, 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::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + + λ_min = v_mean_global[orientation] - c_mean_global + λ_max = v_mean_global[orientation] + c_mean_global + + return λ_min, λ_max +end + +@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector, + equations::LinearizedEulerEquations3D) + @unpack v_mean_global, c_mean_global = equations + + norm_ = norm(normal_direction) + + v_normal = v_mean_global[1] * normal_direction[1] + + v_mean_global[2] * normal_direction[2] + + v_mean_global[3] * normal_direction[3] + + # The v_normals are already scaled by the norm + λ_min = v_normal - c_mean_global * norm_ + λ_max = v_normal + c_mean_global * norm_ + + return λ_min, λ_max +end + +# Convert conservative variables to primitive +@inline cons2prim(u, equations::LinearizedEulerEquations3D) = u +@inline cons2entropy(u, ::LinearizedEulerEquations3D) = u +end # muladd diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index ea7d9193add..8b684b22d09 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -535,6 +535,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_linearizedeuler_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_linearizedeuler_convergence.jl"), + l2=[ + 0.04452389418193219, 0.03688186699434862, + 0.03688186699434861, 0.03688186699434858, + 0.044523894181932186, + ], + linf=[ + 0.2295447498696467, 0.058369658071546704, + 0.05836965807154648, 0.05836965807154648, 0.2295447498696468, + ]) + # 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 # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_tree_3d_linearizedeuler.jl b/test/test_tree_3d_linearizedeuler.jl new file mode 100644 index 00000000000..00f8d62dad9 --- /dev/null +++ b/test/test_tree_3d_linearizedeuler.jl @@ -0,0 +1,35 @@ + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem") + +@testset "Linearized Euler Equations 3D" begin +#! format: noindent + +@trixi_testset "elixir_linearizedeuler_gauss_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gauss_wall.jl"), + l2=[ + 0.020380328336745232, 0.027122442311921492, + 0.02712244231192152, 8.273108096127844e-17, + 0.020380328336745232, + ], + linf=[ + 0.2916021983572774, 0.32763703462270843, + 0.32763703462270855, 1.641012595221666e-15, + 0.2916021983572774, + ], + tspan=(0.0, 1.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 +end diff --git a/test/test_tree_3d_part2.jl b/test/test_tree_3d_part2.jl index 5c7301654ad..4b9da039f98 100644 --- a/test/test_tree_3d_part2.jl +++ b/test/test_tree_3d_part2.jl @@ -22,6 +22,9 @@ isdir(outdir) && rm(outdir, recursive = true) # Compressible Euler with self-gravity include("test_tree_3d_eulergravity.jl") + + # Linearized Euler + include("test_tree_3d_linearizedeuler.jl") end @trixi_testset "Additional tests in 3D" begin diff --git a/test/test_unit.jl b/test/test_unit.jl index 9a4ae2ede13..b23d4aa60e2 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1536,6 +1536,60 @@ end end end +@testset "Equivalent Wave Speed Estimates" begin + @timed_testset "Linearized Euler 3D" begin + equations = LinearizedEulerEquations3D(v_mean_global = (0.42, 0.37, 0.7), + c_mean_global = 1.0, + rho_mean_global = 1.0) + + normal_x = SVector(1.0, 0.0, 0.0) + normal_y = SVector(0.0, 1.0, 0.0) + normal_z = SVector(0.0, 0.0, 1.0) + + u_ll = SVector(0.3, 0.5, -0.7, 0.1, 1.0) + u_rr = SVector(0.5, -0.2, 0.1, 0.2, 5.0) + + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 1, equations), + max_abs_speed_naive(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 2, equations), + max_abs_speed_naive(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 3, equations), + max_abs_speed_naive(u_ll, u_rr, normal_z, + equations))) + + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 1, equations), + min_max_speed_naive(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 2, equations), + min_max_speed_naive(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 3, equations), + min_max_speed_naive(u_ll, u_rr, normal_z, + equations))) + + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 1, equations), + min_max_speed_davis(u_ll, u_rr, normal_x, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 2, equations), + min_max_speed_davis(u_ll, u_rr, normal_y, + equations))) + @test all(isapprox(x, y) + for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 3, equations), + min_max_speed_davis(u_ll, u_rr, normal_z, + equations))) + end +end + @testset "SimpleKronecker" begin N = 3 From 615fd9a44e4aff8cff8ea2bfb8a1e4fb5caec33d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Apr 2024 08:51:41 +0200 Subject: [PATCH 16/17] set version to v0.7.8 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5e1a60b7723..088f936bfda 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.8-pre" +version = "0.7.8" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From dd051566a9079e22dfa792b145c0b50eaafa37fd Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Apr 2024 08:51:55 +0200 Subject: [PATCH 17/17] set development version to v0.7.9-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 088f936bfda..483e9cd8233 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.8" +version = "0.7.9-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"