From 1ba20926010e7882e5ecbc4d45999558a859dea7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Oct 2023 17:11:19 +0200 Subject: [PATCH 01/18] First attempt to enable save solution with time intervals for SimpleIntegratorSSP --- Project.toml | 1 + .../elixir_euler_shockcapturing_subcell.jl | 2 +- src/Trixi.jl | 1 + src/time_integration/methods_SSP.jl | 26 ++++++++++++++----- 4 files changed, 23 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index d318389a6d2..a61928ca33e 100644 --- a/Project.toml +++ b/Project.toml @@ -42,6 +42,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" [weakdeps] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 6b69e4db563..2b9dd466f0d 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -68,7 +68,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 b65d03e7975..67621ad8fce 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -70,6 +70,7 @@ using TriplotBase: TriplotBase using TriplotRecipes: DGTriPseudocolor @reexport using SimpleUnPack: @unpack using SimpleUnPack: @pack! +using DataStructures: BinaryHeap, FasterForward # finite difference SBP operators using SummationByPartsOperators: AbstractDerivativeOperator, diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index a0ed889968a..8fb58213f2b 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -50,17 +50,19 @@ 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()) + push!(tstops_internal, 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 @@ -73,6 +75,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) @@ -84,6 +87,16 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, finalstep::Bool # added for convenience end +""" + add_tstop!(integrator::SimpleIntegratorSSP, t) +Add +""" +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") + push!(integrator.opts.tstops, integrator.tdir * t) +end + # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) function Base.getproperty(integrator::SimpleIntegratorSSP, field::Symbol) if field === :stats @@ -108,8 +121,9 @@ 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) @@ -216,7 +230,7 @@ end # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true - empty!(integrator.opts.tstops) + #empty!(integrator.opts.tstops) end # used for AMR From e05363d74ac684fce05b944029e62872b30939af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Oct 2023 09:44:15 +0200 Subject: [PATCH 02/18] Importing `add_tstop!` --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 67621ad8fce..62d66f62b4d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -37,7 +37,7 @@ 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! using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect From 64e97843e51dc436e889aebcf46f669242f20b1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Oct 2023 11:10:31 +0200 Subject: [PATCH 03/18] First working version of SimpleIntegratorSSP with SaveSolutionCallback using time intervals --- Project.toml | 1 + src/Trixi.jl | 3 ++- src/time_integration/methods_SSP.jl | 16 +++++++++++++--- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index a61928ca33e..28a5b1d66d2 100644 --- a/Project.toml +++ b/Project.toml @@ -43,6 +43,7 @@ Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" [weakdeps] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" diff --git a/src/Trixi.jl b/src/Trixi.jl index 62d66f62b4d..530b61d2ac5 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, add_tstop! + terminate!, remake, add_tstop!, has_tstop, first_tstop +using OrdinaryDiffEq: modify_dt_for_tstops! using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 8fb58213f2b..ac70e3c870b 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -60,7 +60,7 @@ end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward()) - push!(tstops_internal, last(tspan)) + push!(tstops_internal, 2 * last(tspan)) SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback, false, Inf, maxiters, tstops_internal) end @@ -85,6 +85,8 @@ 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 """ @@ -94,9 +96,15 @@ Add 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") + 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 @@ -126,7 +134,7 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; 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)) @@ -169,6 +177,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] @@ -190,7 +200,7 @@ function solve!(integrator::SimpleIntegratorSSP) integrator.iter += 1 integrator.t += integrator.dt - + # handle callbacks if callbacks isa CallbackSet for cb in callbacks.discrete_callbacks From 05fcaaa8d8d7d6f44dfd8a307729b2242d37f3bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Oct 2023 11:16:06 +0200 Subject: [PATCH 04/18] Improved formatting and updated reference solution for test --- src/time_integration/methods_SSP.jl | 13 ++++++++----- test/test_tree_2d_euler.jl | 4 ++-- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index ac70e3c870b..424c25ec398 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -61,8 +61,10 @@ end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward()) push!(tstops_internal, 2 * last(tspan)) - SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback, false, Inf, maxiters, - tstops_internal) + 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 @@ -96,7 +98,7 @@ Add 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") - if length(integrator.opts.tstops) > 1 + if length(integrator.opts.tstops) > 1 pop!(integrator.opts.tstops) end push!(integrator.opts.tstops, integrator.tdir * t) @@ -134,7 +136,8 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p, (prob = ode,), ode.f, alg, SimpleIntegratorSSPOptions(callback, ode.tspan; - kwargs...), false, true, false) + kwargs...), + false, true, false) # resize container resize!(integrator.p, nelements(integrator.p.solver, integrator.p.cache)) @@ -200,7 +203,7 @@ function solve!(integrator::SimpleIntegratorSSP) integrator.iter += 1 integrator.t += integrator.dt - + # handle callbacks if callbacks isa CallbackSet for cb in callbacks.discrete_callbacks diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 1b8a261a60d..e5db2764eaa 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -65,8 +65,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") @trixi_testset "elixir_euler_shockcapturing_subcell.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_subcell.jl"), - l2 = [0.08508147906199143, 0.04510299017724501, 0.045103019801950375, 0.6930704343869766], - linf = [0.31123546471463326, 0.5616274869594462, 0.5619692712224448, 2.88670199345138]) + l2 = [0.08508152653623638, 0.04510301725066843, 0.04510304668512745, 0.6930705064715306], + linf = [0.31136518019691406, 0.5617651935473419, 0.5621200790240503, 2.8866869108596056]) end @trixi_testset "elixir_euler_blast_wave.jl" begin From 71f5c7ad348efe48c10b4a79f505645c0fe9a3ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Oct 2023 12:48:21 +0200 Subject: [PATCH 05/18] Modified initialization of tstops to ensure a stop at the end of the simulation --- src/time_integration/methods_SSP.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 424c25ec398..7865e1d0963 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -60,6 +60,7 @@ end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward()) + push!(tstops_internal, last(tspan)) push!(tstops_internal, 2 * last(tspan)) SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback, false, Inf, From df34188653685d968dee8c8a8651a404a1e2d354 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 23 Oct 2023 11:44:54 +0200 Subject: [PATCH 06/18] Added missing docstring --- src/time_integration/methods_SSP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 7865e1d0963..effe5a85de5 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -94,7 +94,7 @@ end """ add_tstop!(integrator::SimpleIntegratorSSP, t) -Add +Add a time stop during the time integration process. """ function add_tstop!(integrator::SimpleIntegratorSSP, t) integrator.tdir * (t - integrator.t) < zero(integrator.t) && From d04fd98125dee4ff54f544ecdc2c7ce504c8f98d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 23 Oct 2023 12:02:24 +0200 Subject: [PATCH 07/18] Removed OrdinaryDiffEq from Trixi's dependencies --- Project.toml | 1 - src/Trixi.jl | 2 +- src/time_integration/methods_SSP.jl | 24 ++++++++++++++++++++++++ 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 28a5b1d66d2..a61928ca33e 100644 --- a/Project.toml +++ b/Project.toml @@ -43,7 +43,6 @@ Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" [weakdeps] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" diff --git a/src/Trixi.jl b/src/Trixi.jl index 530b61d2ac5..ea96ad694b4 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -38,7 +38,7 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!, AbstractODEIntegrator, init, step!, check_error, get_proposed_dt, set_proposed_dt!, terminate!, remake, add_tstop!, has_tstop, first_tstop -using OrdinaryDiffEq: modify_dt_for_tstops! + using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index effe5a85de5..f13a8b84ac8 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -247,6 +247,30 @@ function terminate!(integrator::SimpleIntegratorSSP) #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 function Base.resize!(integrator::SimpleIntegratorSSP, new_size) resize!(integrator.u, new_size) From 7aa515b5ca79a237067360fb454d1f19defe0041 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 23 Oct 2023 12:24:03 +0200 Subject: [PATCH 08/18] Empty tstops BinaryHeap during the call to `terminate!(integrator::SimpleIntegratorSSP)` --- src/time_integration/methods_SSP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 09d4e91a434..2b598e1ae3b 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -249,7 +249,7 @@ end # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true - #empty!(integrator.opts.tstops) + extract_all!(tstops_internal) end """ From 4322f5b53ec6b56c96114405a8f82553d8b335e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 23 Oct 2023 13:05:12 +0200 Subject: [PATCH 09/18] Fixed bug and added explanatory comments --- src/Trixi.jl | 2 +- src/time_integration/methods_SSP.jl | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index aa947d00510..dd9fda5d324 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -71,7 +71,7 @@ using TriplotBase: TriplotBase using TriplotRecipes: DGTriPseudocolor @reexport using SimpleUnPack: @unpack using SimpleUnPack: @pack! -using DataStructures: BinaryHeap, FasterForward +using DataStructures: BinaryHeap, FasterForward, extract_all! # finite difference SBP operators using SummationByPartsOperators: AbstractDerivativeOperator, diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 2b598e1ae3b..2234293a1ed 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -65,7 +65,10 @@ end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) 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, @@ -99,11 +102,14 @@ end """ add_tstop!(integrator::SimpleIntegratorSSP, t) -Add a time stop during the time integration process. +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 @@ -226,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 @@ -249,7 +259,6 @@ end # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true - extract_all!(tstops_internal) end """ From 772a3889711b82cddaa5ae2e17a851a233e74499 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 23 Oct 2023 15:03:55 +0200 Subject: [PATCH 10/18] Updated Project.toml --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 34ddd3e4a76..2deb14217a2 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.5.47-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" @@ -42,7 +43,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" [weakdeps] Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @@ -53,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" From a3fb5e5ee872bbaf688da75bea9a2559e8eb3e94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 12 Dec 2023 18:44:05 +0100 Subject: [PATCH 11/18] format --- test/test_tree_2d_euler.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 212fded51bc..65899cd5263 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -217,13 +217,13 @@ end 0.08508152653623638, 0.04510301725066843, 0.04510304668512745, - 0.6930705064715306 + 0.6930705064715306, ], linf=[ 0.31136518019691406, 0.5617651935473419, 0.5621200790240503, - 2.8866869108596056 + 2.8866869108596056, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 31d27ad8d6bd3449c374b058f91220de0ef03e64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 19 Dec 2023 15:46:24 +0100 Subject: [PATCH 12/18] Implemented GLM divergence cleaning for multi-ion MHD --- .../tree_2d_dgsem/elixir_mhdmultiion_ec.jl | 9 +- src/callbacks_step/glm_speed_dg.jl | 6 +- src/equations/ideal_mhd_multiion_2d.jl | 108 ++++++++++++++---- 3 files changed, 94 insertions(+), 29 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl index b774abd3047..c2247869c9a 100644 --- a/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl +++ b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl @@ -40,12 +40,17 @@ save_solution = SaveSolutionCallback(dt = 0.1, # interval=100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.5) +cfl = 0.5 + +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) ############################################################################### # run the simulation diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed_dg.jl index 302aae356ab..61426b0d95e 100644 --- a/src/callbacks_step/glm_speed_dg.jl +++ b/src/callbacks_step/glm_speed_dg.jl @@ -7,7 +7,8 @@ function calc_dt_for_cleaning_speed(cfl::Real, mesh, equations::Union{AbstractIdealGlmMhdEquations, - AbstractIdealGlmMhdMulticomponentEquations}, + AbstractIdealGlmMhdMulticomponentEquations, + IdealMhdMultiIonEquations2D}, dg::DG, cache) # compute time step for GLM linear advection equation with c_h=1 for the DG discretization on # Cartesian meshes @@ -20,7 +21,8 @@ end function calc_dt_for_cleaning_speed(cfl::Real, mesh, equations::Union{AbstractIdealGlmMhdEquations, - AbstractIdealGlmMhdMulticomponentEquations}, + AbstractIdealGlmMhdMulticomponentEquations, + IdealMhdMultiIonEquations2D}, dg::DGMulti, cache) rd = dg.basis md = mesh.md diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index ef333d0bb3c..985c260cbe7 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -16,29 +16,31 @@ mutable struct IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real, 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 function IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT, ElectronPressure}(gammas ::SVector{NCOMP, RealT}, charge_to_mass ::SVector{NCOMP, RealT}, electron_pressure - ::ElectronPressure) where + ::ElectronPressure, + c_h::RealT) where {NVARS, NCOMP, RealT <: Real, ElectronPressure} NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value")) - new(gammas, charge_to_mass, electron_pressure) + new(gammas, charge_to_mass, electron_pressure, c_h) end end function IdealMhdMultiIonEquations2D(; gammas, charge_to_mass, - electron_pressure = electron_pressure_zero) + electron_pressure = electron_pressure_zero, + initial_c_h = convert(eltype(gammas), NaN)) _gammas = promote(gammas...) _charge_to_mass = promote(charge_to_mass...) RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass)) - NVARS = length(_gammas) * 5 + 3 + NVARS = length(_gammas) * 5 + 4 NCOMP = length(_gammas) __gammas = SVector(map(RealT, _gammas)) @@ -46,7 +48,8 @@ function IdealMhdMultiIonEquations2D(; gammas, charge_to_mass, return IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT, typeof(electron_pressure)}(__gammas, __charge_to_mass, - electron_pressure) + electron_pressure, + initial_c_h) end @inline function Base.real(::IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT}) where { @@ -66,6 +69,7 @@ function varnames(::typeof(cons2cons), equations::IdealMhdMultiIonEquations2D) tuple("rho_" * string(i), "rho_v1_" * string(i), "rho_v2_" * string(i), "rho_v3_" * string(i), "rho_e_" * string(i))...) end + cons = (cons..., "psi") return cons end @@ -77,6 +81,7 @@ function varnames(::typeof(cons2prim), equations::IdealMhdMultiIonEquations2D) tuple("rho_" * string(i), "v1_" * string(i), "v2_" * string(i), "v3_" * string(i), "p_" * string(i))...) end + prim = (prim..., "psi") return prim end @@ -155,16 +160,18 @@ end # Calculate 1D flux in for a single point @inline function flux(u, orientation::Integer, equations::IdealMhdMultiIonEquations2D) B1, B2, B3 = magnetic_field(u, equations) - + psi = divergence_cleaning_field(u, equations) + v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u, equations) mag_en = 0.5 * (B1^2 + B2^2 + B3^2) + div_clean_energy = 0.5 * psi^2 f = zero(MVector{nvariables(equations), eltype(u)}) if orientation == 1 - f[1] = 0 + f[1] = equations.c_h * psi f[2] = v1_plus * B2 - v2_plus * B1 f[3] = v1_plus * B3 - v3_plus * B1 @@ -176,21 +183,21 @@ end kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - kin_en - mag_en) + p = (gamma - 1) * (rho_e - kin_en - mag_en - div_clean_energy) f1 = rho_v1 f2 = rho_v1 * v1 + p f3 = rho_v1 * v2 f4 = rho_v1 * v3 f5 = (kin_en + gamma * p / (gamma - 1)) * v1 + 2 * mag_en * vk1_plus[k] - - B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + equations.c_h * psi * B1 set_component!(f, k, f1, f2, f3, f4, f5, equations) end - + f[end] = equations.c_h * B1 else #if orientation == 2 f[1] = v2_plus * B1 - v1_plus * B2 - f[2] = 0 + f[2] = equations.c_h * psi f[3] = v2_plus * B3 - v3_plus * B2 for k in eachcomponent(equations) @@ -201,17 +208,18 @@ end kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - kin_en - mag_en) + p = (gamma - 1) * (rho_e - kin_en - mag_en - div_clean_energy) f1 = rho_v2 f2 = rho_v2 * v1 f3 = rho_v2 * v2 + p f4 = rho_v2 * v3 f5 = (kin_en + gamma * p / (gamma - 1)) * v2 + 2 * mag_en * vk2_plus[k] - - B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + equations.c_h * psi * B2 set_component!(f, k, f1, f2, f3, f4, f5, equations) end + f[end] = equations.c_h * B2 end return SVector(f) @@ -271,6 +279,8 @@ The term is composed of three parts # Unpack left and right states to get the magnetic field 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) # Compute important averages B1_avg = 0.5 * (B1_ll + B1_rr) @@ -342,9 +352,14 @@ The term is composed of three parts f4 += charge_ratio_ll[k] * B3_ll * B1_rr f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * B1_rr + # Compute GLM term for the energy + f5 += v1_plus_ll * psi_ll * psi_rr + # Append to the flux vector set_component!(f, k, 0, f2, f3, f4, f5, equations) end + # Compute GLM term for psi + f[end] = v1_plus_ll * psi_rr else #if orientation == 2 # Entries of Powell term for induction equation (already in Trixi's non-conservative form) @@ -386,9 +401,14 @@ The term is composed of three parts f4 += charge_ratio_ll[k] * B3_ll * B2_rr f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * B2_rr + # Compute GLM term for the energy + f5 += v2_plus_ll * psi_ll * psi_rr + # Append to the flux vector set_component!(f, k, 0, f2, f3, f4, f5, equations) end + # Compute GLM term for psi + f[end] = v2_plus_ll * psi_rr end return SVector(f) @@ -407,6 +427,8 @@ The term is composed of three parts # Unpack left and right states to get the magnetic field 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) # Compute important averages mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 @@ -457,11 +479,17 @@ The term is composed of three parts f4 += charge_ratio_ll[k] * B3_ll * B1_rr f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * B1_rr + # Compute GLM term for the energy + f5 += v1_plus_ll * psi_ll * psi_rr + # It's not needed to adjust to Trixi's non-conservative form # Append to the flux vector set_component!(f, k, 0, f2, f3, f4, f5, equations) end + # Compute GLM term for psi + f[end] = v1_plus_ll * psi_rr + else #if orientation == 2 # Entries of Powell term for induction equation (already in Trixi's non-conservative form) f[1] = v1_plus_ll * B2_rr @@ -488,11 +516,16 @@ The term is composed of three parts f4 += charge_ratio_ll[k] * B3_ll * B2_rr f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * B2_rr + # Compute GLM term for the energy + f5 += v2_plus_ll * psi_ll * psi_rr + # It's not needed to adjust to Trixi's non-conservative form # Append to the flux vector set_component!(f, k, 0, f2, f3, f4, f5, equations) end + # Compute GLM term for psi + f[end] = v2_plus_ll * psi_rr end return SVector(f) @@ -515,6 +548,8 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, # Unpack left and right states to get the magnetic field 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) v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll, equations) @@ -533,17 +568,22 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) + psi_avg = 0.5 * (psi_ll + psi_rr) if orientation == 1 + psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) + # Magnetic field components from f^MHD - f6 = 0 + f6 = equations.c_h * psi_avg f7 = v1_plus_avg * B2_avg - v2_plus_avg * B1_avg f8 = v1_plus_avg * B3_avg - v3_plus_avg * B1_avg + f9 = equations.c_h * B1_avg # Start building the flux f[1] = f6 f[2] = f7 f[3] = f8 + f[end] = f9 # Iterate over all components for k in eachcomponent(equations) @@ -563,9 +603,9 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 p_ll = (gammas[k] - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) p_rr = (gammas[k] - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) + (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) beta_ll = 0.5 * rho_ll / p_ll beta_rr = 0.5 * rho_rr / p_rr # for convenience store vk_plus⋅B @@ -614,6 +654,7 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, # MHD part f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v1_plus_mag_avg + B1_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + + f9 * psi_avg - equations.c_h * psi_B1_avg # GLM term + 0.5 * vk1_plus_avg * mag_norm_avg - vk1_plus_avg * B1_avg * B1_avg - vk2_plus_avg * B1_avg * B2_avg - vk3_plus_avg * B1_avg * B3_avg # Additional terms coming from the MHD non-conservative term (momentum eqs) @@ -624,15 +665,19 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, set_component!(f, k, f1, f2, f3, f4, f5, equations) end else #if orientation == 2 + psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) + # Magnetic field components from f^MHD f6 = v2_plus_avg * B1_avg - v1_plus_avg * B2_avg - f7 = 0 + f7 = equations.c_h * psi_avg f8 = v2_plus_avg * B3_avg - v3_plus_avg * B2_avg - + f9 = equations.c_h * B2_avg + # Start building the flux f[1] = f6 f[2] = f7 f[3] = f8 + f[end] = f9 # Iterate over all components for k in eachcomponent(equations) @@ -652,9 +697,9 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 p_ll = (gammas[k] - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) p_rr = (gammas[k] - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) + (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) beta_ll = 0.5 * rho_ll / p_ll beta_rr = 0.5 * rho_rr / p_rr # for convenience store vk_plus⋅B @@ -703,6 +748,7 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, # MHD part f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v2_plus_mag_avg + B2_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + + f9 * psi_avg - equations.c_h * psi_B2_avg # GLM term + 0.5 * vk2_plus_avg * mag_norm_avg - vk1_plus_avg * B2_avg * B1_avg - vk2_plus_avg * B2_avg * B2_avg - vk3_plus_avg * B2_avg * B3_avg # Additional terms coming from the MHD non-conservative term (momentum eqs) @@ -772,6 +818,7 @@ Convert conservative variables to primitive function cons2prim(u, equations::IdealMhdMultiIonEquations2D) @unpack gammas = equations B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(u, equations) prim = zero(MVector{nvariables(equations), eltype(u)}) prim[1] = B1 @@ -786,10 +833,12 @@ function cons2prim(u, equations::IdealMhdMultiIonEquations2D) p = (gammas[k] - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 - + B1 * B1 + B2 * B2 + B3 * B3)) + + B1 * B1 + B2 * B2 + B3 * B3 + + psi * psi)) set_component!(prim, k, rho, v1, v2, v3, p, equations) end + prim[end] = psi return SVector(prim) end @@ -800,6 +849,7 @@ Convert conservative variables to entropy @inline function cons2entropy(u, equations::IdealMhdMultiIonEquations2D) @unpack gammas = equations B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(u, equations) prim = cons2prim(u, equations) entropy = zero(MVector{nvariables(equations), eltype(u)}) @@ -822,6 +872,7 @@ Convert conservative variables to entropy entropy[1] = rho_p_plus * B1 entropy[2] = rho_p_plus * B2 entropy[3] = rho_p_plus * B3 + entropy[end] = rho_p_plus * psi return SVector(entropy) end @@ -832,6 +883,7 @@ Convert primitive to conservative variables @inline function prim2cons(prim, equations::IdealMhdMultiIonEquations2D) @unpack gammas = equations B1, B2, B3 = magnetic_field(prim, equations) + psi = divergence_cleaning_field(prim, equations) cons = zero(MVector{nvariables(equations), eltype(prim)}) cons[1] = B1 @@ -845,10 +897,12 @@ Convert primitive to conservative variables rho_e = p / (gammas[k] - 1.0) + 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * (B1^2 + B2^2 + B3^2) + + 0.5 * psi^2 set_component!(cons, k, rho, rho_v1, rho_v2, rho_v3, rho_e, equations) end + cons[end] = psi return SVector(cons) end @@ -860,6 +914,7 @@ Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoaco @inline function calc_fast_wavespeed(cons, orientation::Integer, equations::IdealMhdMultiIonEquations2D) B1, B2, B3 = magnetic_field(cons, equations) + psi = divergence_cleaning_field(cons, equations) c_f = zero(real(equations)) for k in eachcomponent(equations) @@ -870,7 +925,7 @@ Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoaco v3 = rho_v3 / rho v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2)) + p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) a_square = gamma * p / rho sqrt_rho = sqrt(rho) @@ -953,6 +1008,7 @@ Set the flow variables of component k end magnetic_field(u, equations::IdealMhdMultiIonEquations2D) = SVector(u[1], u[2], u[3]) +divergence_cleaning_field(u, equations::IdealMhdMultiIonEquations2D) = u[end] @inline function density(u, equations::IdealMhdMultiIonEquations2D) rho = zero(real(equations)) @@ -967,6 +1023,8 @@ Computes the sum of the densities times the sum of the pressures """ @inline function density_pressure(u, equations::IdealMhdMultiIonEquations2D) B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(cons, equations) + rho_total = zero(real(equations)) p_total = zero(real(equations)) for k in eachcomponent(equations) @@ -978,7 +1036,7 @@ Computes the sum of the densities times the sum of the pressures v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2)) + p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) rho_total += rho p_total += p From a33ad14b9ae232deda6d118cd5d359d27199c2db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 24 Jan 2024 16:28:14 +0100 Subject: [PATCH 13/18] Added new ES numerical fluxes --- src/Trixi.jl | 1 + src/equations/ideal_mhd_multiion_2d.jl | 268 +++++++++++++++++++++++++ 2 files changed, 269 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index 663c07a29e5..eb65828692f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -178,6 +178,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, + DissipationEntropyStableLump, DissipationEntropyStable, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, min_max_speed_chen_noelle, diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index 985c260cbe7..c508b6a9fe2 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -1043,4 +1043,272 @@ Computes the sum of the densities times the sum of the pressures end return rho_total * p_total end + +""" + DissipationEntropyStableLump(max_abs_speed=max_abs_speed_naive) + +Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed +is estimated as +`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, +defaulting to [`max_abs_speed_naive`](@ref). +""" +struct DissipationEntropyStableLump{MaxAbsSpeed} + max_abs_speed::MaxAbsSpeed +end + +DissipationEntropyStableLump() = DissipationEntropyStableLump(max_abs_speed_naive) + +@inline function (dissipation::DissipationEntropyStableLump)(u_ll, u_rr, + orientation_or_normal_direction, + equations::IdealMhdMultiIonEquations2D) + @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.5 * (B1_ll + B1_rr) + B2_avg = 0.5 * (B2_ll + B2_rr) + B3_avg = 0.5 * (B3_ll + B3_rr) + psi_avg = 0.5 * (psi_ll + psi_rr) + + dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) + + beta_plus_ll = 0 + beta_plus_rr = 0 + # Get the lumped dissipation for all components + 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.5 * rho_ll / p_ll + beta_rr = 0.5 * 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.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) + p_mean = 0.5 * rho_avg / beta_avg + p_star = 0.5 * rho_ln / beta_ln + vel_norm_avg = 0.5 * (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.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) + + h11 = rho_ln + h22 = rho_ln * v1_avg^2 + p_mean + h33 = rho_ln * v2_avg^2 + p_mean + h44 = rho_ln * v3_avg^2 + p_mean + h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln + + vel_avg_norm * p_mean + + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) + + d1 = -0.5 * λ * h11 * (w1_rr - w1_ll) + d2 = -0.5 * λ * h22 * (w2_rr - w2_ll) + d3 = -0.5 * λ * h33 * (w3_rr - w3_ll) + d4 = -0.5 * λ * h44 * (w4_rr - w4_ll) + d5 = -0.5 * λ * 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 + + # Set the magnetic field and psi terms + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) + dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2]) + dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3]) + dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end]) + + return dissipation +end + +function Base.show(io::IO, d::DissipationEntropyStableLump) + print(io, "DissipationEntropyStableLump(", d.max_abs_speed, ")") +end + +""" +DissipationEntropyStable(max_abs_speed=max_abs_speed_naive) + +Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed +is estimated as +`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, +defaulting to [`max_abs_speed_naive`](@ref). +""" +struct DissipationEntropyStable{MaxAbsSpeed} + max_abs_speed::MaxAbsSpeed +end + +DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) + +@inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr, + orientation_or_normal_direction, + equations::IdealMhdMultiIonEquations2D) + @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.5 * (B1_ll + B1_rr) + B2_avg = 0.5 * (B2_ll + B2_rr) + B3_avg = 0.5 * (B3_ll + B3_rr) + psi_avg = 0.5 * (psi_ll + psi_rr) + + dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) + + beta_plus_ll = 0 + beta_plus_rr = 0 + # Get the lumped dissipation for all components + 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.5 * rho_ll / p_ll + beta_rr = 0.5 * 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.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) + p_mean = 0.5 * rho_avg / beta_avg + p_star = 0.5 * rho_ln / beta_ln + vel_norm_avg = 0.5 * (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.5 * 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.5 * λ * (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.5 * λ * (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.5 * λ * (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.5 * λ * (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 + + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) + d5 = -0.5 * λ * (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 + + # Set the magnetic field and psi terms + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + + # diagonal entries + dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) + dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2]) + dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3]) + dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end]) + # Off-diagonal entries + for k in eachcomponent(equations) + _, _, _, _, w5_ll = get_component(k, w_ll, equations) + _, _, _, _, w5_rr = get_component(k, w_rr, equations) + + dissipation[1] -= 0.5 * λ * h_B_psi * B1_avg * (w5_rr - w5_ll) + dissipation[2] -= 0.5 * λ * h_B_psi * B2_avg * (w5_rr - w5_ll) + dissipation[3] -= 0.5 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll) + dissipation[end] -= 0.5 * λ * h_B_psi * psi_avg * (w5_rr - w5_ll) + + ind_E = 3 + (k - 1) * 5 + 5 + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end]) + end + + return dissipation +end + +function Base.show(io::IO, d::DissipationEntropyStable) + print(io, "DissipationEntropyStable(", d.max_abs_speed, ")") +end + end # @muladd From a6b831c85bd516c51ddfcb1e9df1ec663f055500 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 15 Feb 2024 15:37:17 +0100 Subject: [PATCH 14/18] Added missing terms to multi-ion ES flux --- src/equations/ideal_mhd_multiion_2d.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index c508b6a9fe2..2c94b73498f 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -1265,8 +1265,7 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) h53 = h35 h54 = h45 h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln - + vel_avg_norm * p_mean - + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) + + vel_avg_norm * p_mean ) d5 = -0.5 * λ * (h51 * (w1_rr - w1_ll) + h52 * (w2_rr - w2_ll) + h53 * (w3_rr - w3_ll) + @@ -1297,11 +1296,18 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) dissipation[3] -= 0.5 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll) dissipation[end] -= 0.5 * λ * 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 ind_E = 3 + (k - 1) * 5 + 5 dissipation[ind_E] -= 0.5 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1]) dissipation[ind_E] -= 0.5 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2]) dissipation[ind_E] -= 0.5 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3]) dissipation[ind_E] -= 0.5 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end]) + + # Dissipation for the energy equation of all ion species depending on w_5 + for kk in eachcomponent(equations) + ind_E = 3 + (kk - 1) * 5 + 5 + dissipation[ind_E] -= 0.5 * λ * (h_B_psi * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) * (w5_rr - w5_ll) + end end return dissipation From 9fbd286a24e8790d9a1dc8d3ceb26749730a9477 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 15 Feb 2024 15:38:19 +0100 Subject: [PATCH 15/18] Export magnetic_field and divergence_clenaing_field --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index eb65828692f..7f381fd8185 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -220,7 +220,7 @@ export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure, density_product export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, - enstrophy + enstrophy, magnetic_field, divergence_cleaning_field export lake_at_rest_error export ncomponents, eachcomponent export get_component From f0f55a99549e43a28811a3552e120f20b05cac7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 21 Feb 2024 16:10:58 +0100 Subject: [PATCH 16/18] Removed 'lumped' ES flux as it is not really useful --- src/equations/ideal_mhd_multiion_2d.jl | 103 ------------------------- 1 file changed, 103 deletions(-) diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index 2c94b73498f..40d57de7ef3 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -1044,109 +1044,6 @@ Computes the sum of the densities times the sum of the pressures return rho_total * p_total end -""" - DissipationEntropyStableLump(max_abs_speed=max_abs_speed_naive) - -Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed -is estimated as -`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, -defaulting to [`max_abs_speed_naive`](@ref). -""" -struct DissipationEntropyStableLump{MaxAbsSpeed} - max_abs_speed::MaxAbsSpeed -end - -DissipationEntropyStableLump() = DissipationEntropyStableLump(max_abs_speed_naive) - -@inline function (dissipation::DissipationEntropyStableLump)(u_ll, u_rr, - orientation_or_normal_direction, - equations::IdealMhdMultiIonEquations2D) - @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.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - - dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) - - beta_plus_ll = 0 - beta_plus_rr = 0 - # Get the lumped dissipation for all components - 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.5 * rho_ll / p_ll - beta_rr = 0.5 * 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.5 * (rho_ll + rho_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - tau = 1 / (beta_ll + beta_rr) - p_mean = 0.5 * rho_avg / beta_avg - p_star = 0.5 * rho_ln / beta_ln - vel_norm_avg = 0.5 * (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.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) - - h11 = rho_ln - h22 = rho_ln * v1_avg^2 + p_mean - h33 = rho_ln * v2_avg^2 + p_mean - h44 = rho_ln * v3_avg^2 + p_mean - h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln - + vel_avg_norm * p_mean - + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) - - d1 = -0.5 * λ * h11 * (w1_rr - w1_ll) - d2 = -0.5 * λ * h22 * (w2_rr - w2_ll) - d3 = -0.5 * λ * h33 * (w3_rr - w3_ll) - d4 = -0.5 * λ * h44 * (w4_rr - w4_ll) - d5 = -0.5 * λ * 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 - - # Set the magnetic field and psi terms - h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) - dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) - dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2]) - dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3]) - dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end]) - - return dissipation -end - -function Base.show(io::IO, d::DissipationEntropyStableLump) - print(io, "DissipationEntropyStableLump(", d.max_abs_speed, ")") -end - """ DissipationEntropyStable(max_abs_speed=max_abs_speed_naive) From 79ca7c0683b6b35c3c3a56fb6af7eec653c45ac4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 21 Feb 2024 16:12:07 +0100 Subject: [PATCH 17/18] Fixed bug in the discretization of the electron pressure! --- src/equations/ideal_mhd_multiion_2d.jl | 29 +++++++++++++++----------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index 40d57de7ef3..c9e9b64311b 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -291,8 +291,9 @@ The term is composed of three parts mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) # Mean electron pressure - pe_mean = 0.5 * (equations.electron_pressure(u_ll, equations) + - equations.electron_pressure(u_rr, equations)) + pe_ll = equations.electron_pressure(u_ll, equations) + pe_rr = equations.electron_pressure(u_rr, equations) + pe_mean = 0.5 * (pe_ll + pe_rr) # Compute charge ratio of u_ll charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)}) @@ -339,12 +340,13 @@ The term is composed of three parts B3_ll * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) # Adjust non-conservative terms 2 and 3 to Trixi discretization: CHANGE!?! - f2 = 2 * f2 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B1_ll * B1_ll) + f2 = 2 * f2 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B1_ll * B1_ll + pe_ll) f3 = 2 * f3 + charge_ratio_ll[k] * B1_ll * B2_ll f4 = 2 * f4 + charge_ratio_ll[k] * B1_ll * B3_ll - f5 = (2 * f5 - B2_ll * (vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) - - - B3_ll * (vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll)) + f5 = (2 * f5 + - vk1_plus_ll[k] * pe_ll + - B2_ll * (vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) + - B3_ll * (vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll)) # Compute Powell term (already consistent with Trixi's non-conservative discretization) f2 += charge_ratio_ll[k] * B1_ll * B1_rr @@ -389,11 +391,12 @@ The term is composed of three parts # Adjust non-conservative terms 2 and 3 to Trixi discretization: CHANGE!?! f2 = 2 * f2 + charge_ratio_ll[k] * B2_ll * B1_ll - f3 = 2 * f3 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B2_ll * B2_ll) + f3 = 2 * f3 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B2_ll * B2_ll + pe_ll) f4 = 2 * f4 + charge_ratio_ll[k] * B2_ll * B3_ll - f5 = (2 * f5 - B1_ll * (vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) - - - B3_ll * (vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll)) + f5 = (2 * f5 + - vk2_plus_ll[k] * pe_ll + - B1_ll * (vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) + - B3_ll * (vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll)) # Compute Powell term (already consistent with Trixi's non-conservative discretization) f2 += charge_ratio_ll[k] * B1_ll * B2_rr @@ -1047,8 +1050,10 @@ end """ DissipationEntropyStable(max_abs_speed=max_abs_speed_naive) -Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed -is estimated as +Create a local Lax-Friedrichs-type dissipation operator that is provably entropy stable. +See: + * Rueda-Ramírez et al. (2023) +The maximum absolute wave speed is estimated as `max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, defaulting to [`max_abs_speed_naive`](@ref). """ From 9e1dbfff4fcd3b80b2692c39950efa93d76e8a04 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 21 Feb 2024 17:46:37 +0100 Subject: [PATCH 18/18] format --- src/equations/ideal_mhd_multiion_2d.jl | 143 +++++++++++++++---------- 1 file changed, 85 insertions(+), 58 deletions(-) diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index c9e9b64311b..9aeb6f3147c 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -161,7 +161,7 @@ end @inline function flux(u, orientation::Integer, equations::IdealMhdMultiIonEquations2D) B1, B2, B3 = magnetic_field(u, equations) psi = divergence_cleaning_field(u, equations) - + v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u, equations) @@ -190,7 +190,8 @@ end f3 = rho_v1 * v2 f4 = rho_v1 * v3 f5 = (kin_en + gamma * p / (gamma - 1)) * v1 + 2 * mag_en * vk1_plus[k] - - B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + equations.c_h * psi * B1 + B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + + equations.c_h * psi * B1 set_component!(f, k, f1, f2, f3, f4, f5, equations) end @@ -215,7 +216,8 @@ end f3 = rho_v2 * v2 + p f4 = rho_v2 * v3 f5 = (kin_en + gamma * p / (gamma - 1)) * v2 + 2 * mag_en * vk2_plus[k] - - B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + equations.c_h * psi * B2 + B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + + equations.c_h * psi * B2 set_component!(f, k, f1, f2, f3, f4, f5, equations) end @@ -340,13 +342,17 @@ The term is composed of three parts B3_ll * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) # Adjust non-conservative terms 2 and 3 to Trixi discretization: CHANGE!?! - f2 = 2 * f2 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B1_ll * B1_ll + pe_ll) + f2 = 2 * f2 - + charge_ratio_ll[k] * (0.5 * mag_norm_ll - B1_ll * B1_ll + pe_ll) f3 = 2 * f3 + charge_ratio_ll[k] * B1_ll * B2_ll f4 = 2 * f4 + charge_ratio_ll[k] * B1_ll * B3_ll f5 = (2 * f5 - - vk1_plus_ll[k] * pe_ll - - B2_ll * (vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) - - B3_ll * (vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll)) + - + vk1_plus_ll[k] * pe_ll + - + B2_ll * (vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) + - + B3_ll * (vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll)) # Compute Powell term (already consistent with Trixi's non-conservative discretization) f2 += charge_ratio_ll[k] * B1_ll * B1_rr @@ -391,12 +397,16 @@ The term is composed of three parts # Adjust non-conservative terms 2 and 3 to Trixi discretization: CHANGE!?! f2 = 2 * f2 + charge_ratio_ll[k] * B2_ll * B1_ll - f3 = 2 * f3 - charge_ratio_ll[k] * (0.5 * mag_norm_ll - B2_ll * B2_ll + pe_ll) + f3 = 2 * f3 - + charge_ratio_ll[k] * (0.5 * mag_norm_ll - B2_ll * B2_ll + pe_ll) f4 = 2 * f4 + charge_ratio_ll[k] * B2_ll * B3_ll - f5 = (2 * f5 - - vk2_plus_ll[k] * pe_ll - - B1_ll * (vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) - - B3_ll * (vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll)) + f5 = (2 * f5 + - + vk2_plus_ll[k] * pe_ll + - + B1_ll * (vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) + - + B3_ll * (vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll)) # Compute Powell term (already consistent with Trixi's non-conservative discretization) f2 += charge_ratio_ll[k] * B1_ll * B2_rr @@ -606,9 +616,11 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 p_ll = (gammas[k] - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - + 0.5 * psi_ll^2) p_rr = (gammas[k] - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - + 0.5 * psi_rr^2) beta_ll = 0.5 * rho_ll / p_ll beta_rr = 0.5 * rho_rr / p_rr # for convenience store vk_plus⋅B @@ -658,7 +670,8 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v1_plus_mag_avg + B1_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + f9 * psi_avg - equations.c_h * psi_B1_avg # GLM term - + 0.5 * vk1_plus_avg * mag_norm_avg - + + + 0.5 * vk1_plus_avg * mag_norm_avg - vk1_plus_avg * B1_avg * B1_avg - vk2_plus_avg * B1_avg * B2_avg - vk3_plus_avg * B1_avg * B3_avg # Additional terms coming from the MHD non-conservative term (momentum eqs) - @@ -675,7 +688,7 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, f7 = equations.c_h * psi_avg f8 = v2_plus_avg * B3_avg - v3_plus_avg * B2_avg f9 = equations.c_h * B2_avg - + # Start building the flux f[1] = f6 f[2] = f7 @@ -700,9 +713,11 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 p_ll = (gammas[k] - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - + 0.5 * psi_ll^2) p_rr = (gammas[k] - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - + 0.5 * psi_rr^2) beta_ll = 0.5 * rho_ll / p_ll beta_rr = 0.5 * rho_rr / p_rr # for convenience store vk_plus⋅B @@ -752,7 +767,8 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer, f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v2_plus_mag_avg + B2_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + f9 * psi_avg - equations.c_h * psi_B2_avg # GLM term - + 0.5 * vk2_plus_avg * mag_norm_avg - + + + 0.5 * vk2_plus_avg * mag_norm_avg - vk1_plus_avg * B2_avg * B1_avg - vk2_plus_avg * B2_avg * B2_avg - vk3_plus_avg * B2_avg * B3_avg # Additional terms coming from the MHD non-conservative term (momentum eqs) - @@ -928,7 +944,8 @@ Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoaco v3 = rho_v3 / rho v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + p = (gamma - 1) * + (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) a_square = gamma * p / rho sqrt_rho = sqrt(rho) @@ -1039,7 +1056,8 @@ Computes the sum of the densities times the sum of the pressures v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = equations.gammas[k] - p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + p = (gamma - 1) * + (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) rho_total += rho p_total += p @@ -1064,12 +1082,12 @@ end DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) @inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr, - orientation_or_normal_direction, - equations::IdealMhdMultiIonEquations2D) + orientation_or_normal_direction, + equations::IdealMhdMultiIonEquations2D) @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) @@ -1093,7 +1111,7 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) 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) @@ -1111,68 +1129,75 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) v2_avg = 0.5 * (v2_ll + v2_rr) v3_avg = 0.5 * (v3_ll + v3_rr) beta_avg = 0.5 * (beta_ll + beta_rr) - tau = 1 / (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) p_mean = 0.5 * rho_avg / beta_avg p_star = 0.5 * rho_ln / beta_ln vel_norm_avg = 0.5 * (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.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) - + E_bar = p_star / (gammas[k] - 1) + + 0.5 * 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.5 * λ * (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) ) + d1 = -0.5 * λ * + (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.5 * λ * (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) ) + d2 = -0.5 * λ * + (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.5 * λ * (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) ) + d3 = -0.5 * λ * + (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.5 * λ * (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) ) + d4 = -0.5 * λ * + (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.5 * λ * (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) ) + h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln + + + vel_avg_norm * p_mean) + d5 = -0.5 * λ * + (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 @@ -1181,7 +1206,7 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) end # Set the magnetic field and psi terms - h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) # diagonal entries dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) @@ -1208,7 +1233,10 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) # Dissipation for the energy equation of all ion species depending on w_5 for kk in eachcomponent(equations) ind_E = 3 + (kk - 1) * 5 + 5 - dissipation[ind_E] -= 0.5 * λ * (h_B_psi * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) * (w5_rr - w5_ll) + dissipation[ind_E] -= 0.5 * λ * + (h_B_psi * + (B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2)) * + (w5_rr - w5_ll) end end @@ -1218,5 +1246,4 @@ end function Base.show(io::IO, d::DissipationEntropyStable) print(io, "DissipationEntropyStable(", d.max_abs_speed, ")") end - end # @muladd