From 1f31fd47c6c6faf9cba1ab90c6ec7569209e7865 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Dec 2024 09:48:37 -0500 Subject: [PATCH 1/2] Add a Lax-Friedrichs-type entropy-stable dissipation operator and its specialization for multi-ion MHD (#2204) * Added entropy-stable dissipation operator and specialization for multi-ion MHD * Improved commits * Update examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl Co-authored-by: Hendrik Ranocha * Improved comments, reerote test path and reverted real type for multi-ion MHD initial condition to eltype(x) * Update src/equations/ideal_glm_mhd_multiion.jl Co-authored-by: Andrew Winters * Improved comments * Update src/equations/ideal_glm_mhd_multiion.jl Co-authored-by: Andrew Winters * Apply suggestions from code review Co-authored-by: Daniel Doehring * Renamed DissipationEntropyStable to DissipationLaxFriedrichsEntropyVariables * Correct docstring and comment * Update README.md * Update src/equations/ideal_glm_mhd_multiion.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --------- Co-authored-by: Hendrik Ranocha Co-authored-by: Andrew Winters Co-authored-by: Daniel Doehring Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- README.md | 1 + .../tree_2d_dgsem/elixir_mhdmultiion_ec.jl | 7 + src/Trixi.jl | 1 + src/equations/ideal_glm_mhd_multiion.jl | 210 +++++++++++++++++- src/equations/ideal_glm_mhd_multiion_2d.jl | 6 +- src/equations/numerical_fluxes.jl | 40 ++++ test/test_tree_2d_mhdmultiion.jl | 49 +++- test/test_type.jl | 3 + 8 files changed, 312 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 543a0339e41..7397deb4cc4 100644 --- a/README.md +++ b/README.md @@ -54,6 +54,7 @@ installation and postprocessing procedures. Its features include: * Compressible Navier-Stokes equations * Magnetohydrodynamics (MHD) equations * Multi-component compressible Euler and MHD equations + * Multi-ion compressible MHD equations * Linearized Euler and acoustic perturbation equations * Hyperbolic diffusion equations for elliptic problems * Lattice-Boltzmann equations (D2Q9 and D3Q27 schemes) diff --git a/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl index d9b66496b19..491d631f514 100644 --- a/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl +++ b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl @@ -9,8 +9,15 @@ equations = IdealGlmMhdMultiIonEquations2D(gammas = (1.4, 1.667), initial_condition = initial_condition_weak_blast_wave +# Entropy conservative numerical fluxes volume_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal) surface_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal) +# For provably entropy-stable surface fluxes, use +# surface_flux = (FluxPlusDissipation(flux_ruedaramirez_etal, DissipationLaxFriedrichsEntropyVariables()), +# flux_nonconservative_ruedaramirez_etal) +# For a standard local Lax-Friedrichs surface flux, use +# surface_flux = (flux_lax_friedrichs, flux_nonconservative_central) + solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/src/Trixi.jl b/src/Trixi.jl index cb50ca2e7df..f8cc3a1a4d3 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -188,6 +188,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, + DissipationLaxFriedrichsEntropyVariables, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, FluxLMARS, diff --git a/src/equations/ideal_glm_mhd_multiion.jl b/src/equations/ideal_glm_mhd_multiion.jl index c76351a3750..557cb8a12da 100644 --- a/src/equations/ideal_glm_mhd_multiion.jl +++ b/src/equations/ideal_glm_mhd_multiion.jl @@ -9,6 +9,14 @@ have_nonconservative_terms(::AbstractIdealGlmMhdMultiIonEquations) = True() +# Variable names for the multi-ion GLM-MHD equation +# ATTENTION: the variable order for AbstractIdealGlmMhdMultiIonEquations is different than in the reference +# - A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization +# of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. +# [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). +# The first three entries of the state vector `cons[1:3]` are the magnetic field components. After that, we have chunks +# of 5 entries for the hydrodynamic quantities of each ion species. Finally, the last entry `cons[end]` is the divergence +# cleaning field. function varnames(::typeof(cons2cons), equations::AbstractIdealGlmMhdMultiIonEquations) cons = ("B1", "B2", "B3") for i in eachcomponent(equations) @@ -86,7 +94,7 @@ end """ v1, v2, v3, vk1, vk2, vk3 = charge_averaged_velocities(u, - equations::AbstractIdealGlmMhdMultiIonEquations) + equations::AbstractIdealGlmMhdMultiIonEquations) Compute the charge-averaged velocities (`v1`, `v2`, and `v3`) and each ion species' contribution @@ -263,4 +271,204 @@ end return SVector(cons) end + +# Specialization of [`DissipationLaxFriedrichsEntropyVariables`](@ref) for the multi-ion GLM-MHD equations +# For details on the multi-ion entropy Jacobian ``H`` see +# - A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization +# of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. +# [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). +# Since the entropy Jacobian is a sparse matrix, we do not construct it but directly compute the +# action of its product with the jump in the entropy variables. +# +# ATTENTION: the variable order for AbstractIdealGlmMhdMultiIonEquations is different than in the reference above. +# The first three entries of the state vector `u[1:3]` are the magnetic field components. After that, we have chunks +# of 5 entries for the hydrodynamic quantities of each ion species. Finally, the last entry `u[end]` is the divergence +# cleaning field. +@inline function (dissipation::DissipationLaxFriedrichsEntropyVariables)(u_ll, u_rr, + orientation_or_normal_direction, + equations::AbstractIdealGlmMhdMultiIonEquations) + @unpack gammas = equations + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + + w_ll = cons2entropy(u_ll, equations) + w_rr = cons2entropy(u_rr, equations) + prim_ll = cons2prim(u_ll, equations) + prim_rr = cons2prim(u_rr, equations) + B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + # Some global averages + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + + dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) + + beta_plus_ll = 0 + beta_plus_rr = 0 + + # Compute the dissipation for the hydrodynamic quantities of each ion species `k` + ################################################################################# + + # The for loop below fills the entries of `dissipation` that depend on the entries of the diagonal + # blocks ``A_k`` of the entropy Jacobian ``H`` in the given reference (see equations (80)-(82)), + # but the terms that depend on the magnetic field ``B`` and divergence cleaning field ``psi`` are + # excluded here and considered below. In other words, these are the dissipation values that depend + # on the entries of the entropy Jacobian that are marked in blue in Figure 1 of the reference given above. + for k in eachcomponent(equations) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll = get_component(k, prim_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr = get_component(k, prim_rr, equations) + + w1_ll, w2_ll, w3_ll, w4_ll, w5_ll = get_component(k, w_ll, equations) + w1_rr, w2_rr, w3_rr, w4_rr, w5_rr = get_component(k, w_rr, equations) + + # Auxiliary variables + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 + vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 + + # Mean variables + rho_ln = ln_mean(rho_ll, rho_rr) + beta_ln = ln_mean(beta_ll, beta_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + p_star = 0.5f0 * rho_ln / beta_ln + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + vel_avg_norm = v1_avg^2 + v2_avg^2 + v3_avg^2 + E_bar = p_star / (gammas[k] - 1) + + 0.5f0 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) + + h11 = rho_ln + h12 = rho_ln * v1_avg + h13 = rho_ln * v2_avg + h14 = rho_ln * v3_avg + h15 = E_bar + d1 = -0.5f0 * λ * + (h11 * (w1_rr - w1_ll) + + h12 * (w2_rr - w2_ll) + + h13 * (w3_rr - w3_ll) + + h14 * (w4_rr - w4_ll) + + h15 * (w5_rr - w5_ll)) + + h21 = h12 + h22 = rho_ln * v1_avg^2 + p_mean + h23 = h21 * v2_avg + h24 = h21 * v3_avg + h25 = (E_bar + p_mean) * v1_avg + d2 = -0.5f0 * λ * + (h21 * (w1_rr - w1_ll) + + h22 * (w2_rr - w2_ll) + + h23 * (w3_rr - w3_ll) + + h24 * (w4_rr - w4_ll) + + h25 * (w5_rr - w5_ll)) + + h31 = h13 + h32 = h23 + h33 = rho_ln * v2_avg^2 + p_mean + h34 = h31 * v3_avg + h35 = (E_bar + p_mean) * v2_avg + d3 = -0.5f0 * λ * + (h31 * (w1_rr - w1_ll) + + h32 * (w2_rr - w2_ll) + + h33 * (w3_rr - w3_ll) + + h34 * (w4_rr - w4_ll) + + h35 * (w5_rr - w5_ll)) + + h41 = h14 + h42 = h24 + h43 = h34 + h44 = rho_ln * v3_avg^2 + p_mean + h45 = (E_bar + p_mean) * v3_avg + d4 = -0.5f0 * λ * + (h41 * (w1_rr - w1_ll) + + h42 * (w2_rr - w2_ll) + + h43 * (w3_rr - w3_ll) + + h44 * (w4_rr - w4_ll) + + h45 * (w5_rr - w5_ll)) + + h51 = h15 + h52 = h25 + h53 = h35 + h54 = h45 + h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln + + + vel_avg_norm * p_mean) + d5 = -0.5f0 * λ * + (h51 * (w1_rr - w1_ll) + + h52 * (w2_rr - w2_ll) + + h53 * (w3_rr - w3_ll) + + h54 * (w4_rr - w4_ll) + + h55 * (w5_rr - w5_ll)) + + beta_plus_ll += beta_ll + beta_plus_rr += beta_rr + + set_component!(dissipation, k, d1, d2, d3, d4, d5, equations) + end + + # Compute the dissipation related to the magnetic and divergence-cleaning fields + ################################################################################ + + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + + # Dissipation for the magnetic field components due to the diagonal entries of the + # dissipation matrix ``H``. These are the dissipation values that depend on the diagonal + # entries of the entropy Jacobian that are marked in cyan in Figure 1 of the reference given above. + dissipation[1] = -0.5f0 * λ * h_B_psi * (w_rr[1] - w_ll[1]) + dissipation[2] = -0.5f0 * λ * h_B_psi * (w_rr[2] - w_ll[2]) + dissipation[3] = -0.5f0 * λ * h_B_psi * (w_rr[3] - w_ll[3]) + + # Dissipation for the divergence-cleaning field due to the diagonal entry of the + # dissipation matrix ``H``. This dissipation value depends on the single diagonal + # entry of the entropy Jacobian that is marked in red in Figure 1 of the reference given above. + dissipation[end] = -0.5f0 * λ * h_B_psi * (w_rr[end] - w_ll[end]) + + # Dissipation due to the off-diagonal blocks (``B_{off}``) of the dissipation matrix ``H`` and to the entries + # of the block ``A`` that depend on the magnetic field ``B`` and the divergence cleaning field ``psi``. + # See equations (80)-(82) of the given reference. + for k in eachcomponent(equations) + _, _, _, _, w5_ll = get_component(k, w_ll, equations) + _, _, _, _, w5_rr = get_component(k, w_rr, equations) + + # Dissipation for the magnetic field components and divergence cleaning field due to the off-diagonal + # entries of the dissipation matrix ``H`` (block ``B^T`` in equation (80) and Figure 1 of the reference + # given above). + dissipation[1] -= 0.5f0 * λ * h_B_psi * B1_avg * (w5_rr - w5_ll) + dissipation[2] -= 0.5f0 * λ * h_B_psi * B2_avg * (w5_rr - w5_ll) + dissipation[3] -= 0.5f0 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll) + dissipation[end] -= 0.5f0 * λ * h_B_psi * psi_avg * (w5_rr - w5_ll) + + # Dissipation for the energy equation of species `k` depending on `w_1`, `w_2`, `w_3` and `w_end`. These are the + # values of the dissipation that depend on the off-diagonal block ``B`` of the dissipation matrix ``H`` (see equation (80) + # and Figure 1 of the reference given above. + ind_E = 3 + 5 * k # simplified version of 3 + (k - 1) * 5 + 5 + dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1]) + dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2]) + dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3]) + dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end]) + + # Dissipation for the energy equation of all ion species depending on `w_5`. These are the values of the dissipation + # vector that depend on the magnetic and divergence-cleaning field terms of the entries marked with a red cross in + # Figure 1 of the reference given above. + for kk in eachcomponent(equations) + ind_E = 3 + 5 * kk # simplified version of 3 + (kk - 1) * 5 + 5 + dissipation[ind_E] -= 0.5f0 * λ * + (h_B_psi * + (B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2)) * + (w5_rr - w5_ll) + end + end + + return dissipation +end end diff --git a/src/equations/ideal_glm_mhd_multiion_2d.jl b/src/equations/ideal_glm_mhd_multiion_2d.jl index 3bebdefc0e3..fe452665ed9 100644 --- a/src/equations/ideal_glm_mhd_multiion_2d.jl +++ b/src/equations/ideal_glm_mhd_multiion_2d.jl @@ -43,8 +43,8 @@ mutable struct IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real, AbstractIdealGlmMhdMultiIonEquations{2, NVARS, NCOMP} gammas::SVector{NCOMP, RealT} # Heat capacity ratios charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios - electron_pressure::ElectronPressure # Function to compute the electron pressure - c_h::RealT # GLM cleaning speed + electron_pressure::ElectronPressure # Function to compute the electron pressure + c_h::RealT # GLM cleaning speed function IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT, ElectronPressure}(gammas ::SVector{NCOMP, RealT}, @@ -110,7 +110,7 @@ function initial_condition_weak_blast_wave(x, t, # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates - RealT = real(equations) + RealT = eltype(x) inicenter = (0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index ea75b99b7f2..0df67db05a9 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -221,6 +221,46 @@ See [`FluxLaxFriedrichs`](@ref). """ const flux_lax_friedrichs = FluxLaxFriedrichs() +@doc raw""" + DissipationLaxFriedrichsEntropyVariables(max_abs_speed=max_abs_speed_naive) + +Create a local Lax-Friedrichs-type dissipation operator that is provably entropy stable. This operator +must be used together with an entropy-conservative two-point flux function (e.g., `flux_ec`) to yield +an entropy-stable surface flux. The surface flux function can be initialized as: +```julia +flux_es = FluxPlusDissipation(flux_ec, DissipationLaxFriedrichsEntropyVariables()) +``` + +In particular, the numerical flux has the form +```math +f^{\mathrm{ES}} = f^{\mathrm{EC}} - \frac{1}{2} \lambda_{\mathrm{max}} H (w_r - w_l), +``` +where ``f^{\mathrm{EC}}`` is the entropy-conservative two-point flux function (computed with, e.g., `flux_ec`), ``\lambda_{\mathrm{max}}`` +is the maximum wave speed estimated as `max_abs_speed(u_l, u_r, orientation_or_normal_direction, equations)`, +defaulting to [`max_abs_speed_naive`](@ref), ``H`` is a symmetric positive-definite dissipation matrix that +depends on the left and right states `u_l` and `u_r`, and ``(w_r - w_l)`` is the jump in entropy variables. +Ideally, ``H (w_r - w_l) = (u_r - u_l)``, such that the dissipation operator is consistent with the local +Lax-Friedrichs dissipation. + +The entropy-stable dissipation operator is computed with the function +`function (dissipation::DissipationLaxFriedrichsEntropyVariables)(u_l, u_r, orientation_or_normal_direction, equations)`, +which must be specialized for each equation. + +For the derivation of the dissipation matrix for the multi-ion GLM-MHD equations, see: +- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization + of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. + [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). +""" +struct DissipationLaxFriedrichsEntropyVariables{MaxAbsSpeed} + max_abs_speed::MaxAbsSpeed +end + +DissipationLaxFriedrichsEntropyVariables() = DissipationLaxFriedrichsEntropyVariables(max_abs_speed_naive) + +function Base.show(io::IO, d::DissipationLaxFriedrichsEntropyVariables) + print(io, "DissipationLaxFriedrichsEntropyVariables(", d.max_abs_speed, ")") +end + """ FluxHLL(min_max_speed=min_max_speed_davis) diff --git a/test/test_tree_2d_mhdmultiion.jl b/test/test_tree_2d_mhdmultiion.jl index 2ef8b56e9e8..8a9a261c30a 100644 --- a/test/test_tree_2d_mhdmultiion.jl +++ b/test/test_tree_2d_mhdmultiion.jl @@ -6,7 +6,7 @@ using Trixi include("test_trixi.jl") # pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory -EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem") +EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @testset "MHD Multi-ion" begin #! format: noindent @@ -55,6 +55,53 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2 end end +@trixi_testset "Provably entropy-stable LLF-type fluxes for multi-ion GLM-MHD" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"), + l2=[ + 0.017668017558288736, + 0.01779783612885502, + 0.027841673842076285, + 0.015603429086471762, + 0.017849042999817964, + 0.01814196379994667, + 0.005478212889809162, + 0.20585517887094282, + 0.021301245733548135, + 0.03018506565829777, + 0.02938517728342881, + 0.01837279433780041, + 0.11810307914710033, + 0.0002962677911603057 + ], + linf=[ + 0.06594754030722516, + 0.06587779693691242, + 0.09451464686853495, + 0.06787230638663028, + 0.08910065803824378, + 0.08828064474448032, + 0.023647579422062297, + 0.8059383650828509, + 0.1224367642558366, + 0.15930418161523857, + 0.15382860284948224, + 0.08695364286964764, + 0.4949375933243716, + 0.003287251595115295 + ], + surface_flux=(FluxPlusDissipation(flux_ruedaramirez_etal, + DissipationLaxFriedrichsEntropyVariables()), + flux_nonconservative_ruedaramirez_etal)) + # 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_mhdmultiion_ec.jl with local Lax-Friedrichs at the surface" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"), l2=[ diff --git a/test/test_type.jl b/test/test_type.jl index e26b34cdfc6..6f42b7600a4 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1524,6 +1524,7 @@ isdir(outdir) && rm(outdir, recursive = true) one(RealT), one(RealT), one(RealT)) + dissipation_es = DissipationLaxFriedrichsEntropyVariables() orientations = [1, 2] @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == @@ -1547,6 +1548,8 @@ isdir(outdir) && rm(outdir, recursive = true) @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred dissipation_es(u_ll, u_rr, orientation, equations)) == + RealT end @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT From bde0afcf7edcc572c248ee47b9a8f991228aab25 Mon Sep 17 00:00:00 2001 From: Johannes Markert <10619309+jmark@users.noreply.github.com> Date: Wed, 18 Dec 2024 20:26:46 +0100 Subject: [PATCH 2/2] Feature: t8code cubed spherical shell setups (#1918) * Wrap from_abaqus routines. * Implement geometry data transfer from t8code to Trixi. * Updated examples. * Fixed typos. * Applied formatter. * cubed sphere test case, copied from p4est * add baroclinic instability (copy of p4est) * add cubed sphere constructor * Fix indentation. * Fix indentation. * Switching off formatter in two files. * Upgrading T8code.jl. * Fixed examples. * stress different meaning of first argument it refers to level of refinement in lat lon direction, not number of tree as in the p4est version * use lat lon levels * add t8code cubed sphere tests * Remove TODO comments * Relaxing T8code.jl version requirement. * Restricted t8code version requirement. * Removed cubed spherical shell related code. * Removed cubed spherical shell tests. * Including cubed spherical shell setup. * Modified GH ci such that it uses the latest T8code.jl package. * Included cubed spherical shell setups and tests. * Switching order of CI steps. * Changed line in CI. * Increasing code coverage. * Increasing code coverage. * Further increasing code coverage. * Re-introduced tests. * Removed workaround in github workflow. * Applied formatter. * adapt to merged t8code PR * missed calling constructor * Backup. * Add save callback to elixir. * Backup. * Refined code. Make it work in parallel. * Added support for parallelt8codemesh save solution callback. * Applied formatter. * Updated examples and tests. * Applied formatter. * Minor adjustments. * Code refinement. Enabled partitioning after mesh loading. * Applied formatter and fixed typos. * Removed commented out section. * Added missing union type member. * Switching from UInt64 to UInt128 in global interface/mortar id computation. * Switching from UInt64 to UInt128 in global interface/mortar id computation (II). * Adding more tests. * Applied formatter. * Removed Printf. * Incorporated review comments and code polish. * Applied formatter. * Fixed typos. * Update examples/t8code_2d_dgsem/elixir_advection_restart.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update examples/t8code_3d_dgsem/elixir_advection_restart.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/meshes/mesh_io.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update examples/t8code_3d_dgsem/elixir_advection_restart.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Removing last test in t8code 2D MPI to investigate problems in Github CI. * Refactored a bit. * Applied formatter. * Removed commented code. * Added LOG_LEVEL variable. * Added t8code interface simplfication and stitched memory leak. * Applied formatter. * Simplifying finailze behavior for T8codeMesh. * Addeing finalize call to T8codeMesh examples. * Applied formatter. * Update test/test_t8code_3d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update test/test_t8code_3d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * re-add test * remove duplicate * Fixing minor stuff. * bump t8code * Update test/test_t8code_3d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --------- Co-authored-by: Johannes Markert Co-authored-by: Benedict Geihe Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: Hendrik Ranocha Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- Project.toml | 2 +- .../elixir_advection_cubed_sphere.jl | 61 ++++ .../elixir_euler_baroclinic_instability.jl | 299 ++++++++++++++++++ src/meshes/t8code_mesh.jl | 48 ++- src/solvers/dgsem_t8code/containers_3d.jl | 3 +- test/test_t8code_3d.jl | 48 +++ 6 files changed, 457 insertions(+), 4 deletions(-) create mode 100644 examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl create mode 100644 examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl diff --git a/Project.toml b/Project.toml index 296883f5420..95e1da07fd0 100644 --- a/Project.toml +++ b/Project.toml @@ -108,7 +108,7 @@ StaticArrays = "1.5" StrideArrays = "0.1.26" StructArrays = "0.6.11" SummationByPartsOperators = "0.5.41" -T8code = "0.7.2" +T8code = "0.7.4" TimerOutputs = "0.5.7" Triangulate = "2.2" TriplotBase = "0.1" diff --git a/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl b/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl new file mode 100644 index 00000000000..7efffbe52e3 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl @@ -0,0 +1,61 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# 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) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:inside => boundary_condition, + :outside => boundary_condition) + +# Note that the first argument refers to the level of refinement, unlike in for p4est +mesh = Trixi.T8codeMeshCubedSphere(5, 3, 0.5, 0.5; + polydeg = 3, initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +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() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 100, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.2) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + 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/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl new file mode 100644 index 00000000000..d3509069681 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl @@ -0,0 +1,299 @@ +# An idealized baroclinic instability test case +# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. +# +# Note that this elixir can take several hours to run. +# Using 24 threads of an AMD Ryzen Threadripper 3990X (more threads don't speed it up further) +# and `check-bounds=no`, this elixirs takes about one hour to run. +# With 16x16x8 trees per cube face on the same machine, it takes about 28 hours. +# +# References: +# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) +# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores +# https://doi.org/10.1002/qj.2241 + +using OrdinaryDiffEq +using Trixi +using LinearAlgebra + +############################################################################### +# Setup for the baroclinic instability test +gamma = 1.4 +equations = CompressibleEulerEquations3D(gamma) + +# Initial condition for an idealized baroclinic instability test +# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A +function initial_condition_baroclinic_instability(x, t, + equations::CompressibleEulerEquations3D) + lon, lat, r = cartesian_to_sphere(x) + radius_earth = 6.371229e6 + # Make sure that the r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + + # Unperturbed basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Stream function type perturbation + u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) + + u += u_perturbation + v = v_perturbation + + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u - sin(lat) * cos(lon) * v + v2 = cos(lon) * u - sin(lat) * sin(lon) * v + v3 = cos(lat) * v + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +# Steady state for RHS correction below +function steady_state_baroclinic_instability(x, t, equations::CompressibleEulerEquations3D) + lon, lat, r = cartesian_to_sphere(x) + radius_earth = 6.371229e6 + # Make sure that the r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + + # Unperturbed basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u + v2 = cos(lon) * u + v3 = 0.0 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +function cartesian_to_sphere(x) + r = norm(x) + lambda = atan(x[2], x[1]) + if lambda < 0 + lambda += 2 * pi + end + phi = asin(x[3] / r) + + return lambda, phi, r +end + +# Unperturbed balanced steady-state. +# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). +# The other velocity components are zero. +function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + radius_earth = 6.371229e6 # a + half_width_parameter = 2 # b + gravitational_acceleration = 9.80616 # g + k = 3 # k + surface_pressure = 1e5 # p₀ + gas_constant = 287 # R + surface_equatorial_temperature = 310.0 # T₀ᴱ + surface_polar_temperature = 240.0 # T₀ᴾ + lapse_rate = 0.005 # Γ + angular_velocity = 7.29212e-5 # Ω + + # Distance to the center of the Earth + r = z + radius_earth + + # In the paper: T₀ + temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) + # In the paper: A, B, C, H + const_a = 1 / lapse_rate + const_b = (temperature0 - surface_polar_temperature) / + (temperature0 * surface_polar_temperature) + const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / + (surface_equatorial_temperature * surface_polar_temperature) + const_h = gas_constant * temperature0 / gravitational_acceleration + + # In the paper: (r - a) / bH + scaled_z = z / (half_width_parameter * const_h) + + # Temporary variables + temp1 = exp(lapse_rate / temperature0 * z) + temp2 = exp(-scaled_z^2) + + # In the paper: ̃τ₁, ̃τ₂ + tau1 = const_a * lapse_rate / temperature0 * temp1 + + const_b * (1 - 2 * scaled_z^2) * temp2 + tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 + + # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' + inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 + inttau2 = const_c * z * temp2 + + # Temporary variables + temp3 = r / radius_earth * cos(lat) + temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) + + # In the paper: T + temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) + + # In the paper: U, u (zonal wind, first component of spherical velocity) + big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * + (temp3^(k - 1) - temp3^(k + 1)) + temp5 = radius_earth * cos(lat) + u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) + + # Hydrostatic pressure + p = surface_pressure * + exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) + + # Density (via ideal gas law) + rho = p / (gas_constant * temperature) + + return rho, u, p +end + +# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) +function perturbation_stream_function(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + perturbation_radius = 1 / 6 # d₀ / a + perturbed_wind_amplitude = 1.0 # Vₚ + perturbation_lon = pi / 9 # Longitude of perturbation location + perturbation_lat = 2 * pi / 9 # Latitude of perturbation location + pertz = 15000 # Perturbation height cap + + # Great circle distance (d in the paper) divided by a (radius of the Earth) + # because we never actually need d without dividing by a + great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + + cos(perturbation_lat) * cos(lat) * + cos(lon - perturbation_lon)) + + # In the first case, the vertical taper function is per definition zero. + # In the second case, the stream function is per definition zero. + if z > pertz || great_circle_distance_by_a > perturbation_radius + return 0.0, 0.0 + end + + # Vertical tapering of stream function + perttaper = 1.0 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 + + # sin/cos(pi * d / (2 * d_0)) in the paper + sin_, cos_ = sincos(0.5 * pi * great_circle_distance_by_a / perturbation_radius) + + # Common factor for both u and v + factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ + + u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + + cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / + sin(great_circle_distance_by_a) + + v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / + sin(great_circle_distance_by_a) + + return u_perturbation, v_perturbation +end + +@inline function source_terms_baroclinic_instability(u, x, t, + equations::CompressibleEulerEquations3D) + radius_earth = 6.371229e6 # a + gravitational_acceleration = 9.80616 # g + angular_velocity = 7.29212e-5 # Ω + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + r = z + radius_earth + + du1 = zero(eltype(u)) + + # Gravity term + temp = -gravitational_acceleration * radius_earth^2 / r^3 + du2 = temp * u[1] * x[1] + du3 = temp * u[1] * x[2] + du4 = temp * u[1] * x[3] + du5 = temp * (u[2] * x[1] + u[3] * x[2] + u[4] * x[3]) + + # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4] + du2 -= -2 * angular_velocity * u[3] + du3 -= 2 * angular_velocity * u[2] + + return SVector(du1, du2, du3, du4, du5) +end + +############################################################################### +# Start of the actual elixir, semidiscretization of the problem + +initial_condition = initial_condition_baroclinic_instability + +boundary_conditions = Dict(:inside => boundary_condition_slip_wall, + :outside => boundary_condition_slip_wall) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340) +volume_flux = flux_kennedy_gruber +solver = DGSEM(polydeg = 5, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# For optimal results, use (16, 8) here +trees_per_cube_face = (8, 4) +mesh = Trixi.T8codeMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0, + polydeg = 5, initial_refinement_level = 0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_baroclinic_instability, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10 * 24 * 60 * 60.0) # time in seconds for 10 days + +# Save RHS of the steady state and subtract it in every RHS evaluation. +# This trick preserves the steady state exactly (to machine rounding errors, of course). +# Otherwise, this elixir produces entirely unusable results for a resolution of 8x8x4 cells +# per cube face with a polydeg of 3. +# With this trick, even the polydeg 3 simulation produces usable (although badly resolved) results, +# and most of the grid imprinting in higher polydeg simulation is eliminated. +# +# See https://github.com/trixi-framework/Trixi.jl/issues/980 for more information. +u_steady_state = compute_coefficients(steady_state_baroclinic_instability, tspan[1], semi) +# Use a `let` block for performance (otherwise du_steady_state will be a global variable) +let du_steady_state = similar(u_steady_state) + # Save RHS of the steady state + Trixi.rhs!(du_steady_state, u_steady_state, semi, tspan[1]) + + global function corrected_rhs!(du, u, semi, t) + # Normal RHS evaluation + Trixi.rhs!(du, u, semi, t) + # Correct by subtracting the steady-state RHS + Trixi.@trixi_timeit Trixi.timer() "rhs correction" begin + # Use Trixi.@threaded for threaded performance + Trixi.@threaded for i in eachindex(du) + du[i] -= du_steady_state[i] + end + end + end +end +u0 = compute_coefficients(tspan[1], semi) +ode = ODEProblem(corrected_rhs!, u0, tspan, semi) + +summary_callback = SummaryCallback() + +analysis_interval = 5000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 5000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution) + +############################################################################### +# run the simulation + +# Use a Runge-Kutta method with automatic (error based) time step size control +# Enable threading of the RK method for better performance on multiple threads +sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1.0e-6, + reltol = 1.0e-6, + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index fa5b1d8f81a..4c85c543e16 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -598,7 +598,7 @@ struct AbaqusFile{NDIMS} end """ - T8codeMesh(meshfile::String, NDIMS; kwargs...) + T8codeMesh(filepath::String, NDIMS; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from either a Gmsh mesh file (`.msh`) or Abaqus mesh file (`.inp`) which is determined @@ -783,6 +783,52 @@ function t8_cmesh_new_from_connectivity(connectivity::Ptr{p8est_connectivity}, c return t8_cmesh_new_from_p8est(connectivity, comm, 0) end +""" +T8codeMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; + polydeg, RealT=Float64, initial_refinement_level=0) + +Construct a cubed spherical shell of given inner radius and thickness as `T8codeMesh` with +`6 * trees_per_face_dimension^2 * layers` trees. The mesh will have two boundaries, +`:inside` and `:outside`. + +# Arguments +- `lat_lon_levels_per_face_dimension::Integer`: number of trees per patch in longitudinal + and latitudinal direction given as level of + refinement. +- `layers::Integer`: the number of trees in the third local dimension of each face, i.e., + the number of layers of the shell. +- `inner_radius::Float64`: Radius of the inner side of the shell. +- `thickness::Float64`: Thickness of the shell. The outer radius will be + `inner_radius + thickness`. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the + simulation starts. +""" +function T8codeMeshCubedSphere(lat_lon_levels_per_face_dimension, layers, inner_radius, + thickness; + polydeg, RealT = Float64, initial_refinement_level = 0) + NDIMS = 3 + cmesh = t8_cmesh_new_cubed_spherical_shell(inner_radius, thickness, + lat_lon_levels_per_face_dimension, + layers, mpi_comm()) + do_face_ghost = mpi_isparallel() + scheme = t8_scheme_new_default_cxx() + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, + mpi_comm()) + + num_trees = t8_cmesh_get_num_trees(cmesh) + boundary_names = fill(Symbol("---"), 2 * NDIMS, num_trees) + for itree in 1:num_trees + boundary_names[5, itree] = :inside + boundary_names[6, itree] = :outside + end + + return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = polydeg) +end + struct adapt_callback_passthrough adapt_callback::Function user_data::Any diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl index 981dabc397f..a4eb86ebad3 100644 --- a/src/solvers/dgsem_t8code/containers_3d.jl +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -64,8 +64,7 @@ function calc_node_coordinates!(node_coordinates, current_index += 1), matrix1, matrix2, matrix3, view(mesh.tree_node_coordinates, :, :, :, :, - global_itree + 1), - tmp1) + global_itree + 1), tmp1) end end diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index e07cba1f576..82bec769f28 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -111,6 +111,22 @@ mkdir(outdir) end end + # This test differs from the one in `test_p4est_3d.jl` in the latitudinal and + # longitudinal dimensions. + @trixi_testset "elixir_advection_cubed_sphere.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"), + l2=[0.002006918015656413], + linf=[0.027655117058380085]) + # 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 + # This test is identical to the one in `test_p4est_3d.jl`. @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), @@ -336,6 +352,38 @@ mkdir(outdir) end end + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_baroclinic_instability.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_baroclinic_instability.jl"), + l2=[ + 6.725093801700048e-7, + 0.00021710076010951073, + 0.0004386796338203878, + 0.00020836270267103122, + 0.07601887903440395 + ], + linf=[ + 1.9107530539574924e-5, + 0.02980358831035801, + 0.048476331898047564, + 0.02200137344113612, + 4.848310144356219 + ], + tspan=(0.0, 1e2), + # Decrease tolerance of adaptive time stepping to get similar results across different systems + abstol=1.0e-9, reltol=1.0e-9, + coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI + # 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_euler_weak_blast_wave_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), l2=[