diff --git a/Project.toml b/Project.toml index 539dafc303..a834ddeea7 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.6.5-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -52,6 +53,7 @@ TrixiMakieExt = "Makie" [compat] CodeTracking = "1.0.5" ConstructionBase = "1.3" +DataStructures = "0.18.15" DiffEqCallbacks = "2.25" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" @@ -84,7 +86,7 @@ StaticArrays = "1" StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" -T8code = "0.4.3" +T8code = "0.4.3, 0.5" TimerOutputs = "0.5.7" Triangulate = "2.0" TriplotBase = "0.1" diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 282805a0e0..44e63a0872 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -67,7 +67,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 100, +save_solution = SaveSolutionCallback(dt = 0.1, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) diff --git a/src/Trixi.jl b/src/Trixi.jl index e7b849e264..e18b2f6415 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -37,7 +37,8 @@ using SciMLBase: CallbackSet, DiscreteCallback, import SciMLBase: get_du, get_tmp_cache, u_modified!, AbstractODEIntegrator, init, step!, check_error, get_proposed_dt, set_proposed_dt!, - terminate!, remake + terminate!, remake, add_tstop!, has_tstop, first_tstop + using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect @@ -70,6 +71,7 @@ using TriplotBase: TriplotBase using TriplotRecipes: DGTriPseudocolor @reexport using SimpleUnPack: @unpack using SimpleUnPack: @pack! +using DataStructures: BinaryHeap, FasterForward, extract_all! # finite difference SBP operators using SummationByPartsOperators: AbstractDerivativeOperator, diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 8ddb0dcd08..6338e04c3e 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -461,6 +461,7 @@ end # Convert conservative variables to entropy @inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations1D) @unpack cv, gammas, gas_constants = equations + rho_v1, rho_e = u rho = density(u, equations) @@ -480,21 +481,86 @@ end s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p T = (rho_e - 0.5 * rho * v_square) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant * - ((gamma - s) / - (gamma - 1.0) - - (0.5 * v_square * - rho_p)) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 2])) - + v1^2 / (2 * T)) for i in eachcomponent(equations)) w1 = gas_constant * v1 * rho_p - w2 = gas_constant * (-1.0 * rho_p) + w2 = gas_constant * (-rho_p) entrop_other = SVector{2, real(equations)}(w1, w2) return vcat(entrop_other, entrop_rho) end +# Convert entropy variables to conservative variables +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations1D) + @unpack gammas, gas_constants, cv, cp = equations + T = -1 / w[2] + v1 = w[1] * T + cons_rho = SVector{ncomponents(equations), real(equations)}(exp(1 / + gas_constants[i] * + (-cv[i] * + log(-w[2]) - + cp[i] + w[i + 2] - + 0.5 * w[1]^2 / + w[2])) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + gamma = help1 / help2 + u2 = p / (gamma - 1) + 0.5 * rho * v1^2 + cons_other = SVector{2, real(equations)}(u1, u2) + return vcat(cons_other, cons_rho) +end + +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations1D) + @unpack cv, gammas, gas_constants = equations + rho_v1, rho_e = u + rho = density(u, equations) + T = temperature(u, equations) + + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2])) + end + + return total_entropy +end + +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations1D) + @unpack cv, gammas, gas_constants = equations + + rho_v1, rho_e = u + + rho = density(u, equations) + help1 = zero(rho) + + for i in eachcomponent(equations) + help1 += u[i + 2] * cv[i] + end + + v1 = rho_v1 / rho + v_square = v1^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + + return T +end + """ totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 940d88b1aa..60fce222f2 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -665,22 +665,57 @@ end s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p T = (rho_e - 0.5 * rho * v_square) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant * - ((gamma - s) / - (gamma - 1.0) - - (0.5 * v_square * - rho_p)) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) for i in eachcomponent(equations)) w1 = gas_constant * v1 * rho_p w2 = gas_constant * v2 * rho_p - w3 = gas_constant * rho_p * (-1) + w3 = gas_constant * (-rho_p) entrop_other = SVector{3, real(equations)}(w1, w2, w3) return vcat(entrop_other, entrop_rho) end +# Convert entropy variables to conservative variables +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) + @unpack gammas, gas_constants, cp, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end + # Convert primitive to conservative variables @inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D) @unpack cv, gammas = equations @@ -700,6 +735,39 @@ end return vcat(cons_other, cons_rho) end +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + rho = density(u, equations) + T = temperature(u, equations) + + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) + end + + return total_entropy +end + +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + help1 = zero(rho) + + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + + return T +end + """ totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index d2f22679c4..f3b09b01e9 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -204,6 +204,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator2N, dt) integrator.dt = dt end +# Required e.g. for `glm_speed_callback` +function get_proposed_dt(integrator::SimpleIntegrator2N) + return integrator.dt +end + # stop the time integration function terminate!(integrator::SimpleIntegrator2N) integrator.finalstep = true diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index b0ce593051..7b70466606 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/src/time_integration/methods_3Sstar.jl @@ -282,6 +282,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator3Sstar, dt) integrator.dt = dt end +# Required e.g. for `glm_speed_callback` +function get_proposed_dt(integrator::SimpleIntegrator3Sstar) + return integrator.dt +end + # stop the time integration function terminate!(integrator::SimpleIntegrator3Sstar) integrator.finalstep = true diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index dbb9e51121..9d1e06488b 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -55,17 +55,25 @@ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 -mutable struct SimpleIntegratorSSPOptions{Callback} +mutable struct SimpleIntegratorSSPOptions{Callback, TStops} callback::Callback # callbacks; used in Trixi adaptive::Bool # whether the algorithm is adaptive; ignored dtmax::Float64 # ignored maxiters::Int # maximal number of time steps - tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored + tstops::TStops # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) - SimpleIntegratorSSPOptions{typeof(callback)}(callback, false, Inf, maxiters, - [last(tspan)]) + tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward()) + # We add last(tspan) to make sure that the time integration stops at the end time + push!(tstops_internal, last(tspan)) + # We add 2 * last(tspan) because add_tstop!(integrator, t) is only called by DiffEqCallbacks.jl if tstops contains a time that is larger than t + # (https://github.com/SciML/DiffEqCallbacks.jl/blob/025dfe99029bd0f30a2e027582744528eb92cd24/src/iterative_and_periodic.jl#L92) + push!(tstops_internal, 2 * last(tspan)) + SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback, + false, Inf, + maxiters, + tstops_internal) end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77 @@ -78,6 +86,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, du::uType r0::uType t::RealT + tdir::RealT dt::RealT # current time step dtcache::RealT # ignored iter::Int # current number of time steps (iteration) @@ -87,8 +96,29 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, alg::Alg opts::SimpleIntegratorSSPOptions finalstep::Bool # added for convenience + dtchangeable::Bool + force_stepfail::Bool end +""" + add_tstop!(integrator::SimpleIntegratorSSP, t) +Add a time stop during the time integration process. +This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution. +""" +function add_tstop!(integrator::SimpleIntegratorSSP, t) + integrator.tdir * (t - integrator.t) < zero(integrator.t) && + error("Tried to add a tstop that is behind the current time. This is strictly forbidden") + # We need to remove the first entry of tstops when a new entry is added. + # Otherwise, the simulation gets stuck at the previous tstop and dt is adjusted to zero. + if length(integrator.opts.tstops) > 1 + pop!(integrator.opts.tstops) + end + push!(integrator.opts.tstops, integrator.tdir * t) +end + +has_tstop(integrator::SimpleIntegratorSSP) = !isempty(integrator.opts.tstops) +first_tstop(integrator::SimpleIntegratorSSP) = first(integrator.opts.tstops) + # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) function Base.getproperty(integrator::SimpleIntegratorSSP, field::Symbol) if field === :stats @@ -113,11 +143,13 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; du = similar(u) r0 = similar(u) t = first(ode.tspan) + tdir = sign(ode.tspan[end] - ode.tspan[1]) iter = 0 - integrator = SimpleIntegratorSSP(u, du, r0, t, dt, zero(dt), iter, ode.p, + integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p, (prob = ode,), ode.f, alg, SimpleIntegratorSSPOptions(callback, ode.tspan; - kwargs...), false) + kwargs...), + false, true, false) # resize container resize!(integrator.p, nelements(integrator.p.solver, integrator.p.cache)) @@ -160,6 +192,8 @@ function solve!(integrator::SimpleIntegratorSSP) terminate!(integrator) end + modify_dt_for_tstops!(integrator) + @. integrator.r0 = integrator.u for stage in eachindex(alg.c) t_stage = integrator.t + integrator.dt * alg.c[stage] @@ -198,6 +232,10 @@ function solve!(integrator::SimpleIntegratorSSP) end end + # Empty the tstops array. + # This cannot be done in terminate!(integrator::SimpleIntegratorSSP) because DiffEqCallbacks.PeriodicCallbackAffect would return at error. + extract_all!(integrator.opts.tstops) + for stage_callback in alg.stage_callbacks finalize_callback(stage_callback, integrator.p) end @@ -226,7 +264,30 @@ end # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true - empty!(integrator.opts.tstops) +end + +""" + modify_dt_for_tstops!(integrator::SimpleIntegratorSSP) +Modify the time-step size to match the time stops specified in integrator.opts.tstops. +To avoid adding OrdinaryDiffEq to Trixi's dependencies, this routine is a copy of +https://github.com/SciML/OrdinaryDiffEq.jl/blob/d76335281c540ee5a6d1bd8bb634713e004f62ee/src/integrators/integrator_utils.jl#L38-L54 +""" +function modify_dt_for_tstops!(integrator::SimpleIntegratorSSP) + if has_tstop(integrator) + tdir_t = integrator.tdir * integrator.t + tdir_tstop = first_tstop(integrator) + if integrator.opts.adaptive + integrator.dt = integrator.tdir * + min(abs(integrator.dt), abs(tdir_tstop - tdir_t)) # step! to the end + elseif iszero(integrator.dtcache) && integrator.dtchangeable + integrator.dt = integrator.tdir * abs(tdir_tstop - tdir_t) + elseif integrator.dtchangeable && !integrator.force_stepfail + # always try to step! with dtcache, but lower if a tstop + # however, if force_stepfail then don't set to dtcache, and no tstop worry + integrator.dt = integrator.tdir * + min(abs(integrator.dtcache), abs(tdir_tstop - tdir_t)) # step! to the end + end + end end # used for AMR diff --git a/test/test_tree_1d_eulermulti.jl b/test/test_tree_1d_eulermulti.jl index bd86de928e..b6c79ce03d 100644 --- a/test/test_tree_1d_eulermulti.jl +++ b/test/test_tree_1d_eulermulti.jl @@ -2,136 +2,153 @@ module TestExamples1DEulerMulti using Test using Trixi +using ForwardDiff include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Compressible Euler Multicomponent" begin -#! format: noindent + @trixi_testset "Testing entropy2cons and cons2entropy" begin + using ForwardDiff + gammas = (1.3272378792562836, 1.5269959187969864, 1.8362285750521512, + 1.0409061360276926, 1.4652015053812224, 1.3626493264184423) + gas_constants = (1.817636851910076, 6.760820475922636, 5.588953939749113, + 6.31574782981543, 3.362932038038397, 3.212779569399733) + equations = CompressibleEulerMulticomponentEquations1D(gammas = SVector{length(gammas)}(gammas...), + gas_constants = SVector{length(gas_constants)}(gas_constants...)) + u = [-1.4632513788889214, 0.9908786980927811, 0.2909066990257628, + 0.6256623915420473, 0.4905882754313441, 0.14481800501749112, + 1.0333532872771651, 0.6805599818745411] + w = cons2entropy(u, equations) + # test that the entropy variables match the gradients of the total entropy + @test w ≈ ForwardDiff.gradient(u -> Trixi.total_entropy(u, equations), u) + # test that `entropy2cons` is the inverse of `cons2entropy` + @test entropy2cons(w, equations) ≈ u + end -@trixi_testset "elixir_eulermulti_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2=[0.15330089521538684, 0.4417674632047301, - 0.016888510510282385, 0.03377702102056477, - 0.06755404204112954], - linf=[0.29130548795961864, 0.8847009003152357, - 0.034686525099975274, 0.06937305019995055, - 0.1387461003999011]) - # 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 + @trixi_testset "elixir_eulermulti_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), + l2=[0.15330089521538684, 0.4417674632047301, + 0.016888510510282385, 0.03377702102056477, + 0.06755404204112954], + linf=[0.29130548795961864, 0.8847009003152357, + 0.034686525099975274, 0.06937305019995055, + 0.1387461003999011]) + # 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 -@trixi_testset "elixir_eulermulti_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2=[ - 0.1522380497572071, - 0.43830846465313206, - 0.03907262116499431, - 0.07814524232998862, - ], - linf=[ - 0.24939193075537294, - 0.7139395740052739, - 0.06324208768391237, - 0.12648417536782475, - ]) - # 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 + @trixi_testset "elixir_eulermulti_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), + l2=[ + 0.1522380497572071, + 0.43830846465313206, + 0.03907262116499431, + 0.07814524232998862, + ], + linf=[ + 0.24939193075537294, + 0.7139395740052739, + 0.06324208768391237, + 0.12648417536782475, + ]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2=[ - 8.575236038539227e-5, - 0.00016387804318585358, - 1.9412699303977585e-5, - 3.882539860795517e-5, - ], - linf=[ - 0.00030593277277124464, - 0.0006244803933350696, - 7.253121435135679e-5, - 0.00014506242870271358, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 8.575236038539227e-5, + 0.00016387804318585358, + 1.9412699303977585e-5, + 3.882539860795517e-5, + ], + linf=[ + 0.00030593277277124464, + 0.0006244803933350696, + 7.253121435135679e-5, + 0.00014506242870271358, + ]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[1.8983933794407234e-5, 6.207744299844731e-5, - 1.5466205761868047e-6, 3.0932411523736094e-6, - 6.186482304747219e-6, 1.2372964609494437e-5], - linf=[0.00012014372605895218, 0.0003313207215800418, - 6.50836791016296e-6, 1.301673582032592e-5, - 2.603347164065184e-5, 5.206694328130368e-5]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[1.8983933794407234e-5, 6.207744299844731e-5, + 1.5466205761868047e-6, 3.0932411523736094e-6, + 6.186482304747219e-6, 1.2372964609494437e-5], + linf=[0.00012014372605895218, 0.0003313207215800418, + 6.50836791016296e-6, 1.301673582032592e-5, + 2.603347164065184e-5, 5.206694328130368e-5]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[1.888450477353845e-5, 5.4910600482795386e-5, - 9.426737161533622e-7, 1.8853474323067245e-6, - 3.770694864613449e-6, 7.541389729226898e-6], - linf=[0.00011622351152063004, 0.0003079221967086099, - 3.2177423254231563e-6, 6.435484650846313e-6, - 1.2870969301692625e-5, 2.574193860338525e-5], - volume_flux=flux_chandrashekar) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[1.888450477353845e-5, 5.4910600482795386e-5, + 9.426737161533622e-7, 1.8853474323067245e-6, + 3.770694864613449e-6, 7.541389729226898e-6], + linf=[0.00011622351152063004, 0.0003079221967086099, + 3.2177423254231563e-6, 6.435484650846313e-6, + 1.2870969301692625e-5, 2.574193860338525e-5], + volume_flux=flux_chandrashekar) + # 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 -@trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_two_interacting_blast_waves.jl"), - l2=[1.288867611915533, 82.71335258388848, 0.00350680272313187, - 0.013698784353152794, - 0.019179518517518084], - linf=[29.6413044707026, 1322.5844802186496, 0.09191919374782143, - 0.31092970966717925, - 0.4417989757182038], - tspan=(0.0, 0.0001)) - # 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 + @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_two_interacting_blast_waves.jl"), + l2=[1.288867611915533, 82.71335258388848, 0.00350680272313187, + 0.013698784353152794, + 0.019179518517518084], + linf=[29.6413044707026, 1322.5844802186496, 0.09191919374782143, + 0.31092970966717925, + 0.4417989757182038], + tspan=(0.0, 0.0001)) + # 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 -end end # module diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 04a295537a..65899cd526 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -214,16 +214,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_subcell.jl"), l2=[ - 0.08508147906199143, - 0.04510299017724501, - 0.045103019801950375, - 0.6930704343869766, + 0.08508152653623638, + 0.04510301725066843, + 0.04510304668512745, + 0.6930705064715306, ], linf=[ - 0.31123546471463326, - 0.5616274869594462, - 0.5619692712224448, - 2.88670199345138, + 0.31136518019691406, + 0.5617651935473419, + 0.5621200790240503, + 2.8866869108596056, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 30d52b37b9..7c4a4e722e 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -8,234 +8,253 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") @testset "Compressible Euler Multicomponent" begin -#! format: noindent + @trixi_testset "Testing entropy2cons and cons2entropy" begin + using ForwardDiff + gammas = (1.1546412974182538, 1.1171560258914812, 1.097107661471476, + 1.0587601652669245, 1.6209889683979308, 1.6732209755396386, + 1.2954303574165822) + gas_constants = (5.969461071171914, 3.6660802003290183, 6.639008614675539, + 8.116604827140456, 6.190706056680031, 1.6795013743693712, + 2.197737590916966) + equations = CompressibleEulerMulticomponentEquations2D(gammas = SVector{length(gammas)}(gammas...), + gas_constants = SVector{length(gas_constants)}(gas_constants...)) + u = [-1.7433292819144075, 0.8844413258376495, 0.6050737175812364, + 0.8261998359817043, 1.0801186290896465, 0.505654488367698, + 0.6364415555805734, 0.851669392285058, 0.31219606420306223, + 1.0930477805612038] + w = cons2entropy(u, equations) + # test that the entropy variables match the gradients of the total entropy + @test w ≈ ForwardDiff.gradient(u -> Trixi.total_entropy(u, equations), u) + # test that `entropy2cons` is the inverse of `cons2entropy` + @test entropy2cons(w, equations) ≈ u + end -# NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the -# simulations are set up with dimensional states. For example, the reference pressure in SI -# units is 101325 Pa, i.e., pressure has values of O(10^5) + # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the + # simulations are set up with dimensional states. For example, the reference pressure in SI + # units is 101325 Pa, i.e., pressure has values of O(10^5) -@trixi_testset "elixir_eulermulti_shock_bubble.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"), - l2=[ - 73.78467629094177, - 0.9174752929795251, - 57942.83587826468, - 0.1828847253029943, - 0.011127037850925347, - ], - linf=[ - 196.81051991521073, - 7.8456811648529605, - 158891.88930113698, - 0.811379581519794, - 0.08011973559187913, - ], - tspan=(0.0, 0.001)) - # 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 + @trixi_testset "elixir_eulermulti_shock_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"), + l2=[ + 73.78467629094177, + 0.9174752929795251, + 57942.83587826468, + 0.1828847253029943, + 0.011127037850925347, + ], + linf=[ + 196.81051991521073, + 7.8456811648529605, + 158891.88930113698, + 0.811379581519794, + 0.08011973559187913, + ], + tspan=(0.0, 0.001)) + # 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 -@trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin - rm("out/deviations.txt", force = true) - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), - l2=[ - 81.52845664909304, - 2.5455678559421346, - 63229.190712645846, - 0.19929478404550321, - 0.011068604228443425, - ], - linf=[ - 249.21708417382013, - 40.33299887640794, - 174205.0118831558, - 0.6881458768113586, - 0.11274401158173972, - ], - initial_refinement_level=3, - tspan=(0.0, 0.001), - output_directory="out") - lines = readlines("out/deviations.txt") - @test lines[1] == "# iter, simu_time, rho1_min, rho2_min" - @test startswith(lines[end], "1") - # 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)) < 15000 + @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin + rm("out/deviations.txt", force = true) + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), + l2=[ + 81.52845664909304, + 2.5455678559421346, + 63229.190712645846, + 0.19929478404550321, + 0.011068604228443425, + ], + linf=[ + 249.21708417382013, + 40.33299887640794, + 174205.0118831558, + 0.6881458768113586, + 0.11274401158173972, + ], + initial_refinement_level=3, + tspan=(0.0, 0.001), + output_directory="out") + lines = readlines("out/deviations.txt") + @test lines[1] == "# iter, simu_time, rho1_min, rho2_min" + @test startswith(lines[end], "1") + # 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)) < 15000 + end end -end -@trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"), - l2=[ - 73.10832638093902, - 1.4599215762968585, - 57176.014861335476, - 0.17812843581838675, - 0.010123079422717837, - ], - linf=[ - 214.50568817511956, - 25.40392579616452, - 152862.41011222568, - 0.564195553101797, - 0.0956331651771212, - ], - initial_refinement_level=3, - tspan=(0.0, 0.001)) - # 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)) < 15000 + @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"), + l2=[ + 73.10832638093902, + 1.4599215762968585, + 57176.014861335476, + 0.17812843581838675, + 0.010123079422717837, + ], + linf=[ + 214.50568817511956, + 25.40392579616452, + 152862.41011222568, + 0.564195553101797, + 0.0956331651771212, + ], + initial_refinement_level=3, + tspan=(0.0, 0.001)) + # 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)) < 15000 + end end -end -@trixi_testset "elixir_eulermulti_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2=[ - 0.050182236154087095, - 0.050189894464434635, - 0.2258715597305131, - 0.06175171559771687, - ], - linf=[ - 0.3108124923284472, - 0.3107380389947733, - 1.054035804988521, - 0.29347582879608936, - ]) - # 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 + @trixi_testset "elixir_eulermulti_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), + l2=[ + 0.050182236154087095, + 0.050189894464434635, + 0.2258715597305131, + 0.06175171559771687, + ], + linf=[ + 0.3108124923284472, + 0.3107380389947733, + 1.054035804988521, + 0.29347582879608936, + ]) + # 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 -@trixi_testset "elixir_eulermulti_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2=[ - 0.0496546258404055, - 0.04965550099933263, - 0.22425206549856372, - 0.004087155041747821, - 0.008174310083495642, - 0.016348620166991283, - 0.032697240333982566, - ], - linf=[ - 0.2488251110766228, - 0.24832493304479406, - 0.9310354690058298, - 0.017452870465607374, - 0.03490574093121475, - 0.0698114818624295, - 0.139622963724859, - ]) - # 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 + @trixi_testset "elixir_eulermulti_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), + l2=[ + 0.0496546258404055, + 0.04965550099933263, + 0.22425206549856372, + 0.004087155041747821, + 0.008174310083495642, + 0.016348620166991283, + 0.032697240333982566, + ], + linf=[ + 0.2488251110766228, + 0.24832493304479406, + 0.9310354690058298, + 0.017452870465607374, + 0.03490574093121475, + 0.0698114818624295, + 0.139622963724859, + ]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2=[ - 0.00012290225488326508, - 0.00012290225488321876, - 0.00018867397906337653, - 4.8542321753649044e-5, - 9.708464350729809e-5, - ], - linf=[ - 0.0006722819239133315, - 0.0006722819239128874, - 0.0012662292789555885, - 0.0002843844182700561, - 0.0005687688365401122, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 0.00012290225488326508, + 0.00012290225488321876, + 0.00018867397906337653, + 4.8542321753649044e-5, + 9.708464350729809e-5, + ], + linf=[ + 0.0006722819239133315, + 0.0006722819239128874, + 0.0012662292789555885, + 0.0002843844182700561, + 0.0005687688365401122, + ]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[ - 2.2661773867001696e-6, - 2.266177386666318e-6, - 6.593514692980009e-6, - 8.836308667348217e-7, - 1.7672617334696433e-6, - ], - linf=[ - 1.4713170997993075e-5, - 1.4713170997104896e-5, - 5.115618808515521e-5, - 5.3639516094383666e-6, - 1.0727903218876733e-5, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[ + 2.2661773867001696e-6, + 2.266177386666318e-6, + 6.593514692980009e-6, + 8.836308667348217e-7, + 1.7672617334696433e-6, + ], + linf=[ + 1.4713170997993075e-5, + 1.4713170997104896e-5, + 5.115618808515521e-5, + 5.3639516094383666e-6, + 1.0727903218876733e-5, + ]) + # 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 -@trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[ - 1.8621737639352465e-6, - 1.862173764098385e-6, - 5.942585713809631e-6, - 6.216263279534722e-7, - 1.2432526559069443e-6, - ], - linf=[ - 1.6235495582606063e-5, - 1.6235495576388814e-5, - 5.854523678827661e-5, - 5.790274858807898e-6, - 1.1580549717615796e-5, - ], - volume_flux=flux_chandrashekar) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[ + 1.8621737639352465e-6, + 1.862173764098385e-6, + 5.942585713809631e-6, + 6.216263279534722e-7, + 1.2432526559069443e-6, + ], + linf=[ + 1.6235495582606063e-5, + 1.6235495576388814e-5, + 5.854523678827661e-5, + 5.790274858807898e-6, + 1.1580549717615796e-5, + ], + volume_flux=flux_chandrashekar) + # 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 -end end # module