From a9b3caf7066c67bc4489b487c642e72fa1846671 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 08:50:16 +0100 Subject: [PATCH 01/21] initial fix of the c_h computation a la Fluxo. Still ugly and involves logic in the stepsize callback --- src/callbacks_step/stepsize_dg2d.jl | 22 +++++++++++++++------- src/equations/ideal_glm_mhd_2d.jl | 5 ----- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 93f097aa2d..1b5b3cdecc 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -5,13 +5,6 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) - # FIXME: This should be implemented properly using another callback - # or something else, cf. - # https://github.com/trixi-framework/Trixi.jl/issues/257 - if equations isa AbstractIdealGlmMhdEquations - equations.c_h = zero(equations.c_h) - end - for element in eachelement(dg, cache) max_λ1 = max_λ2 = zero(max_scaled_speed) for j in eachnode(dg), i in eachnode(dg) @@ -24,6 +17,21 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2)) end + if equations isa AbstractIdealGlmMhdEquations + # 1) compute time step for ONLY the GLM linear advection equation with c_h = 1 + # 2) must be redone each time due to the AMR possibility + # OBS! possibly could use the max_dt function below somehow for the constant wave speed case + max_scaled_speed_for_c_h = nextfloat(zero(t)) + for element in eachelement(dg, cache) + inv_jacobian = cache.elements.inverse_jacobian[element] + max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) + end + c_h_deltat = 2 / (nnodes(dg) * max_scaled_speed_for_c_h) + # c_h is proportional to its own time step divided by the complete mhd time step + # OBS! scaled by 1/2 for safety + equations.c_h = 0.5 * c_h_deltat / (2 / (nnodes(dg) * max_scaled_speed)) + end + return 2 / (nnodes(dg) * max_scaled_speed) end diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 15e22f3ffd..d4461f7db2 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -546,11 +546,6 @@ end cf_y_direction = calc_fast_wavespeed(u, 2, equations) cf_max = max(cf_x_direction, cf_y_direction) - # FIXME: This should be implemented properly using another callback - # or something else, cf. - # https://github.com/trixi-framework/Trixi.jl/issues/257 - equations.c_h = max(equations.c_h, cf_max) # GLM cleaning speed = c_f - return abs(v1) + cf_x_direction, abs(v2) + cf_y_direction end From d509ef2e68d6a9d31433bc9c29a038a8dd5ae2a2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 09:12:29 +0100 Subject: [PATCH 02/21] updated comments with known issue about c_h time sclaing. Tests will fail need updated --- src/callbacks_step/stepsize_dg2d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 1b5b3cdecc..0434406360 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -21,6 +21,7 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, # 1) compute time step for ONLY the GLM linear advection equation with c_h = 1 # 2) must be redone each time due to the AMR possibility # OBS! possibly could use the max_dt function below somehow for the constant wave speed case + # OBS! This needs to be scaled by the CFL number for consistency max_scaled_speed_for_c_h = nextfloat(zero(t)) for element in eachelement(dg, cache) inv_jacobian = cache.elements.inverse_jacobian[element] From 114ba563e6ac2188afafb8c13408c548b1511cf1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 10:55:40 +0100 Subject: [PATCH 03/21] made glm speed update into a new callback. still needs cleanup and test updating --- examples/2d/elixir_mhd_alfven_wave.jl | 11 ++-- src/Trixi.jl | 2 +- src/auxiliary/precompile.jl | 2 +- src/callbacks_step/callbacks_step.jl | 3 ++ src/callbacks_step/glm_speed.jl | 72 +++++++++++++++++++++++++++ src/callbacks_step/stepsize_dg2d.jl | 30 +++++------ 6 files changed, 99 insertions(+), 21 deletions(-) create mode 100644 src/callbacks_step/glm_speed.jl diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl index dc1ed0e738..504e7de8be 100644 --- a/examples/2d/elixir_mhd_alfven_wave.jl +++ b/examples/2d/elixir_mhd_alfven_wave.jl @@ -36,7 +36,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_anal extra_analysis_integrals=(entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity)) - + alive_callback = AliveCallback(analysis_interval=analysis_interval) save_solution = SaveSolutionCallback(interval=10, @@ -46,11 +46,14 @@ save_solution = SaveSolutionCallback(interval=10, stepsize_callback = StepsizeCallback(cfl=1.0) -callbacks = CallbackSet(summary_callback, - analysis_callback, +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) diff --git a/src/Trixi.jl b/src/Trixi.jl index 1923ddf71c..a4fa9dc3b9 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -109,7 +109,7 @@ export SemidiscretizationEulerGravity, ParametersEulerGravity, export SummaryCallback, SteadyStateCallback, AMRCallback, StepsizeCallback, SaveRestartCallback, SaveSolutionCallback, AnalysisCallback, AliveCallback, - TrivialCallback + TrivialCallback, GlmSpeedCallback export load_mesh, load_time diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index c95e1d2030..0e8902e8de 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -230,6 +230,7 @@ function _precompile_manual_() @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:analysis_interval,),Tuple{Int}},Type{AliveCallback}}) for RealT in (Float64,) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:cfl,),Tuple{RealT}},Type{StepsizeCallback}}) + @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:glm_scale,),Tuple{RealT}},Type{GlmSpeedCallback}}) end @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_final_restart),Tuple{Int,Bool}},Type{SaveRestartCallback}}) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_initial_solution, :save_final_solution, :solution_variables),Tuple{Int,Bool,Bool,Symbol}},Type{SaveSolutionCallback}}) @@ -273,4 +274,3 @@ function _precompile_manual_() return nothing end - diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index d64519407e..9879fded93 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -36,6 +36,7 @@ end # in a time step loop (however, callbacks are always executed *after* a step, thus it comes near # the end here) # * `StepsizeCallback` must come after AMR to accomodate potential changes in the minimum cell size +# * `GlmSpeedCallback` must come after computing time step size because it affects the value of c_h include("summary.jl") include("steady_state.jl") include("analysis.jl") @@ -44,6 +45,8 @@ include("save_restart.jl") include("save_solution.jl") include("amr.jl") include("stepsize.jl") +include("glm_speed.jl") + # The `TrivialCallback` purposely does nothing: It allows to quickly disable specific callbacks # when using `trixi_include` or `test_trixi_include` diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl new file mode 100644 index 0000000000..7cc9bab2a2 --- /dev/null +++ b/src/callbacks_step/glm_speed.jl @@ -0,0 +1,72 @@ + +""" + GlmSpeedCallback(; glm_scale=0.5) + +Update the divergence cleaning wave speed c_h according to the time step +computed in StepsizeCallback for the ideal GLM-MHD equations. +""" +mutable struct GlmSpeedCallback{RealT} + glm_scale::RealT +end + + +function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:GlmSpeedCallback} + glm_speed_callback = cb.affect! + @unpack glm_scale = glm_speed_callback + print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:GlmSpeedCallback} + if get(io, :compact, false) + show(io, cb) + else + glm_speed_callback = cb.affect! + + setup = [ + "GLM wave speed scaling" => glm_speed_callback.glm_scale, + ] + summary_box(io, "GlmSpeedCallback", setup) + end +end + + +function GlmSpeedCallback(; glm_scale::Real=0.5) + # when is the callback activated + condition = (u, t, integrator) -> true + + glm_speed_callback = GlmSpeedCallback(glm_scale) + + DiscreteCallback(condition, glm_speed_callback, + save_positions=(false,false), + initialize=initialize!) +end + + +function initialize!(cb::DiscreteCallback{Condition,Affect!}, u, t, integrator) where {Condition, Affect!<:GlmSpeedCallback} + cb.affect!(integrator) +end + + +# This method is called as callback after the StepsizeCallback during the time integration. +@inline function (glm_speed_callback::GlmSpeedCallback)(integrator) + + dt = integrator.dtcache + semi = integrator.p + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + @unpack glm_scale = glm_speed_callback + + # compute time step for ONLY the GLM linear advection equation with c_h = 1 + # must be redone each time due to the AMR possibility + max_scaled_speed_for_c_h = nextfloat(zero(dt)) + for element in eachelement(solver, cache) + inv_jacobian = cache.elements.inverse_jacobian[element] + max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) + end + cfl = 1.0 # FIXME: This is a hack beacuse I was not sure how to get access to this information from StepsizeCallback + c_h_deltat = cfl * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) + + # c_h is proportional to its own time step divided by the complete mhd time step + equations.c_h = 0.5 * c_h_deltat / dt + + return nothing +end diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 0434406360..f72d034ad8 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -17,21 +17,21 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2)) end - if equations isa AbstractIdealGlmMhdEquations - # 1) compute time step for ONLY the GLM linear advection equation with c_h = 1 - # 2) must be redone each time due to the AMR possibility - # OBS! possibly could use the max_dt function below somehow for the constant wave speed case - # OBS! This needs to be scaled by the CFL number for consistency - max_scaled_speed_for_c_h = nextfloat(zero(t)) - for element in eachelement(dg, cache) - inv_jacobian = cache.elements.inverse_jacobian[element] - max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) - end - c_h_deltat = 2 / (nnodes(dg) * max_scaled_speed_for_c_h) - # c_h is proportional to its own time step divided by the complete mhd time step - # OBS! scaled by 1/2 for safety - equations.c_h = 0.5 * c_h_deltat / (2 / (nnodes(dg) * max_scaled_speed)) - end + # if equations isa AbstractIdealGlmMhdEquations + # # 1) compute time step for ONLY the GLM linear advection equation with c_h = 1 + # # 2) must be redone each time due to the AMR possibility + # # OBS! possibly could use the max_dt function below somehow for the constant wave speed case + # # OBS! This needs to be scaled by the CFL number for consistency + # max_scaled_speed_for_c_h = nextfloat(zero(t)) + # for element in eachelement(dg, cache) + # inv_jacobian = cache.elements.inverse_jacobian[element] + # max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) + # end + # c_h_deltat = 2 / (nnodes(dg) * max_scaled_speed_for_c_h) + # # c_h is proportional to its own time step divided by the complete mhd time step + # # OBS! scaled by 1/2 for safety + # equations.c_h = 0.5 * c_h_deltat / (2 / (nnodes(dg) * max_scaled_speed)) + # end return 2 / (nnodes(dg) * max_scaled_speed) end From e440803a791f7aaaab15da59118eba291352a143 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 11:32:07 +0100 Subject: [PATCH 04/21] now callback is dimension agnostic --- src/callbacks_step/glm_speed.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 7cc9bab2a2..03d819ea47 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -60,7 +60,7 @@ end max_scaled_speed_for_c_h = nextfloat(zero(dt)) for element in eachelement(solver, cache) inv_jacobian = cache.elements.inverse_jacobian[element] - max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) + max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, ndims(semi.equations) * inv_jacobian) end cfl = 1.0 # FIXME: This is a hack beacuse I was not sure how to get access to this information from StepsizeCallback c_h_deltat = cfl * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) From 291d413d0d38d0059742d34f7e066f8870b816d0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 11:45:52 +0100 Subject: [PATCH 05/21] glm callback now stores the cfl number for convenience. fixed bug where glm_scale was unused --- examples/2d/elixir_mhd_alfven_wave.jl | 5 +++-- src/auxiliary/precompile.jl | 2 +- src/callbacks_step/glm_speed.jl | 15 ++++++++------- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl index 504e7de8be..f5be6ffb6c 100644 --- a/examples/2d/elixir_mhd_alfven_wave.jl +++ b/examples/2d/elixir_mhd_alfven_wave.jl @@ -44,9 +44,10 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.0) +cfl_number = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl_number) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 0e8902e8de..d5b6f7c501 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -230,7 +230,7 @@ function _precompile_manual_() @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:analysis_interval,),Tuple{Int}},Type{AliveCallback}}) for RealT in (Float64,) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:cfl,),Tuple{RealT}},Type{StepsizeCallback}}) - @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:glm_scale,),Tuple{RealT}},Type{GlmSpeedCallback}}) + @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:glm_scale, :cfl_scale),Tuple{RealT,RealT}},Type{GlmSpeedCallback}}) end @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_final_restart),Tuple{Int,Bool}},Type{SaveRestartCallback}}) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_initial_solution, :save_final_solution, :solution_variables),Tuple{Int,Bool,Bool,Symbol}},Type{SaveSolutionCallback}}) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 03d819ea47..1ac423956c 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -5,7 +5,8 @@ Update the divergence cleaning wave speed c_h according to the time step computed in StepsizeCallback for the ideal GLM-MHD equations. """ -mutable struct GlmSpeedCallback{RealT} +mutable struct GlmSpeedCallback{RealT<:Real} + cfl_scale::RealT glm_scale::RealT end @@ -16,6 +17,7 @@ function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condi print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ")") end + function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:GlmSpeedCallback} if get(io, :compact, false) show(io, cb) @@ -30,11 +32,11 @@ function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Af end -function GlmSpeedCallback(; glm_scale::Real=0.5) +function GlmSpeedCallback(; glm_scale=0.5, cfl_scale=1.0) # when is the callback activated condition = (u, t, integrator) -> true - glm_speed_callback = GlmSpeedCallback(glm_scale) + glm_speed_callback = GlmSpeedCallback(glm_scale, cfl_scale) DiscreteCallback(condition, glm_speed_callback, save_positions=(false,false), @@ -53,7 +55,7 @@ end dt = integrator.dtcache semi = integrator.p mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - @unpack glm_scale = glm_speed_callback + @unpack glm_scale, cfl_scale = glm_speed_callback # compute time step for ONLY the GLM linear advection equation with c_h = 1 # must be redone each time due to the AMR possibility @@ -62,11 +64,10 @@ end inv_jacobian = cache.elements.inverse_jacobian[element] max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, ndims(semi.equations) * inv_jacobian) end - cfl = 1.0 # FIXME: This is a hack beacuse I was not sure how to get access to this information from StepsizeCallback - c_h_deltat = cfl * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) + c_h_deltat = cfl_scale * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) # c_h is proportional to its own time step divided by the complete mhd time step - equations.c_h = 0.5 * c_h_deltat / dt + equations.c_h = glm_scale * c_h_deltat / dt return nothing end From 615fefa7aba9625daeac3964db779aed1648bcaa Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 11:48:00 +0100 Subject: [PATCH 06/21] removed commented code from stepsize callback --- src/callbacks_step/stepsize_dg2d.jl | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index f72d034ad8..56c453a88c 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -17,22 +17,6 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2)) end - # if equations isa AbstractIdealGlmMhdEquations - # # 1) compute time step for ONLY the GLM linear advection equation with c_h = 1 - # # 2) must be redone each time due to the AMR possibility - # # OBS! possibly could use the max_dt function below somehow for the constant wave speed case - # # OBS! This needs to be scaled by the CFL number for consistency - # max_scaled_speed_for_c_h = nextfloat(zero(t)) - # for element in eachelement(dg, cache) - # inv_jacobian = cache.elements.inverse_jacobian[element] - # max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, 2 * inv_jacobian) - # end - # c_h_deltat = 2 / (nnodes(dg) * max_scaled_speed_for_c_h) - # # c_h is proportional to its own time step divided by the complete mhd time step - # # OBS! scaled by 1/2 for safety - # equations.c_h = 0.5 * c_h_deltat / (2 / (nnodes(dg) * max_scaled_speed)) - # end - return 2 / (nnodes(dg) * max_scaled_speed) end From 1c6dd879dcbf0153dedf7721c2124ef753cb84f8 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 12:12:23 +0100 Subject: [PATCH 07/21] fixed bug in the order of arguments within the struct --- src/callbacks_step/glm_speed.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 1ac423956c..ec1d2934f8 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -1,13 +1,13 @@ """ - GlmSpeedCallback(; glm_scale=0.5) + GlmSpeedCallback(; glm_scale=0.5, cfl_scale=1.0) Update the divergence cleaning wave speed c_h according to the time step computed in StepsizeCallback for the ideal GLM-MHD equations. """ mutable struct GlmSpeedCallback{RealT<:Real} - cfl_scale::RealT glm_scale::RealT + cfl_scale::RealT end @@ -57,8 +57,7 @@ end mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack glm_scale, cfl_scale = glm_speed_callback - # compute time step for ONLY the GLM linear advection equation with c_h = 1 - # must be redone each time due to the AMR possibility + # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) max_scaled_speed_for_c_h = nextfloat(zero(dt)) for element in eachelement(solver, cache) inv_jacobian = cache.elements.inverse_jacobian[element] @@ -66,7 +65,7 @@ end end c_h_deltat = cfl_scale * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) - # c_h is proportional to its own time step divided by the complete mhd time step + # c_h is proportional to its own time step divided by the complete MHD time step equations.c_h = glm_scale * c_h_deltat / dt return nothing From 41a7562360a6005979267afdb8a1e56334eed505 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 12:30:42 +0100 Subject: [PATCH 08/21] added the glm callback to all the MHD examples --- examples/2d/elixir_mhd_alfven_wave.jl | 1 - examples/2d/elixir_mhd_alfven_wave_mortar.jl | 12 ++++++++---- examples/2d/elixir_mhd_blast_wave.jl | 14 +++++++++----- examples/2d/elixir_mhd_ec.jl | 12 ++++++++---- examples/2d/elixir_mhd_orszag_tang.jl | 14 +++++++++----- examples/2d/elixir_mhd_rotor.jl | 14 +++++++++----- 6 files changed, 43 insertions(+), 24 deletions(-) diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl index f5be6ffb6c..c164608296 100644 --- a/examples/2d/elixir_mhd_alfven_wave.jl +++ b/examples/2d/elixir_mhd_alfven_wave.jl @@ -57,7 +57,6 @@ callbacks = CallbackSet(summary_callback, glm_speed_callback) - ############################################################################### # run the simulation diff --git a/examples/2d/elixir_mhd_alfven_wave_mortar.jl b/examples/2d/elixir_mhd_alfven_wave_mortar.jl index bb001e4f52..4ae9a9597c 100644 --- a/examples/2d/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/2d/elixir_mhd_alfven_wave_mortar.jl @@ -47,13 +47,17 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.0) +cfl_number = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl_number) + +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) callbacks = CallbackSet(summary_callback, - analysis_callback, - alive_callback, + analysis_callback, + alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) ############################################################################### # run the simulation diff --git a/examples/2d/elixir_mhd_blast_wave.jl b/examples/2d/elixir_mhd_blast_wave.jl index 1041760c0d..2f6f58a923 100644 --- a/examples/2d/elixir_mhd_blast_wave.jl +++ b/examples/2d/elixir_mhd_blast_wave.jl @@ -63,14 +63,18 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -stepsize_callback = StepsizeCallback(cfl=0.3) +cfl_number = 0.3 +stepsize_callback = StepsizeCallback(cfl=cfl_number) -callbacks = CallbackSet(summary_callback, - analysis_callback, +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - amr_callback, - stepsize_callback) + amr_callback, + stepsize_callback, + glm_speed_callback) ############################################################################### # run the simulation diff --git a/examples/2d/elixir_mhd_ec.jl b/examples/2d/elixir_mhd_ec.jl index 106541add8..7a0dff6a65 100644 --- a/examples/2d/elixir_mhd_ec.jl +++ b/examples/2d/elixir_mhd_ec.jl @@ -41,13 +41,17 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.0) +cfl_number = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl_number) -callbacks = CallbackSet(summary_callback, - analysis_callback, +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) ############################################################################### diff --git a/examples/2d/elixir_mhd_orszag_tang.jl b/examples/2d/elixir_mhd_orszag_tang.jl index 178755868d..4746c3c4dc 100644 --- a/examples/2d/elixir_mhd_orszag_tang.jl +++ b/examples/2d/elixir_mhd_orszag_tang.jl @@ -63,14 +63,18 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -stepsize_callback = StepsizeCallback(cfl=0.4) +cfl_number = 0.4 +stepsize_callback = StepsizeCallback(cfl=cfl_number) -callbacks = CallbackSet(summary_callback, - analysis_callback, +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - amr_callback, - stepsize_callback) + amr_callback, + stepsize_callback, + glm_speed_callback) ############################################################################### # run the simulation diff --git a/examples/2d/elixir_mhd_rotor.jl b/examples/2d/elixir_mhd_rotor.jl index 35dbf0319a..173b76300f 100644 --- a/examples/2d/elixir_mhd_rotor.jl +++ b/examples/2d/elixir_mhd_rotor.jl @@ -64,14 +64,18 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -stepsize_callback = StepsizeCallback(cfl=0.35) +cfl_number = 0.35 +stepsize_callback = StepsizeCallback(cfl=cfl_number) -callbacks = CallbackSet(summary_callback, - analysis_callback, +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - amr_callback, - stepsize_callback) + amr_callback, + stepsize_callback, + glm_speed_callback) ############################################################################### # run the simulation From 97a71eeb1358b495eaa4ae06b780b350ace61beb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 12:57:29 +0100 Subject: [PATCH 09/21] all tests updated and pass locally --- test/test_examples_2d.jl | 12 ++++++------ test/test_examples_2d_mhd.jl | 20 ++++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_examples_2d.jl b/test/test_examples_2d.jl index db8a0813bf..a92f9e2026 100644 --- a/test/test_examples_2d.jl +++ b/test/test_examples_2d.jl @@ -133,23 +133,23 @@ end # GLM-MHD @testset "elixir_mhd_alfven_wave.jl one step with initial_condition_constant" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2 = [7.144325530681224e-17, 2.107107684242054e-16, 5.025035342841136e-16, 3.9063728231270855e-17, 8.448122434419674e-15, 3.9171737639099993e-16, 2.404057005740705e-16, 3.6588423152083e-17, 1.609894978924592e-16], - linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 6.384151957784351e-16], + l2 = [7.144325530681224e-17, 2.123397983547417e-16, 5.061138912500049e-16, 3.6588423152083e-17, 8.449816179702522e-15, 3.9171737639099993e-16, 2.445565690318772e-16, 3.6588423152083e-17, 9.971153407737885e-17], + linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 4.779017148551244e-16], maxiters = 1, initial_condition = initial_condition_constant) end @testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2 = [1.2430371080579117, 1.7996415471431686, 1.6901513554005372, 0.0, 2.263340343580966, 0.21272304014984506, 0.23334102161201117, 0.0, 0.002328993566253826], - linf = [10.402314125114891, 14.060210214680058, 15.559978821022856, 0.0, 16.72422832990639, 1.3007692563104454, 1.412328778460149, 0.0, 0.05937649008864308], + l2 = [1.2429643014570362, 1.7996363685163963, 1.6899889382208215, 0.0, 2.263014495582255, 0.2127948919559293, 0.23341167012961367, 0.0, 0.003341061230142962], + linf = [10.401662016997985, 14.061042946011094, 15.556382737749237, 0.0, 16.714999099689578, 1.3250449624793987, 1.4148467737020096, 0.0, 0.0878998114279951], tspan = (0.0, 0.05)) end @testset "elixir_mhd_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"), - l2 = [0.17570357318216484, 3.853646257410666, 2.47457711714592, 0.0, 355.29948057316955, 2.3505873160544257, 1.3965840971478956, 0.0, 0.02331813093276944], - linf = [1.5824262876996142, 44.191703566463, 12.87587680235509, 0.0, 2236.7827386920094, 13.005033422940064, 9.01906588206558, 0.0, 0.3569041393901721], + l2 = [0.17569823349539196, 3.853292514951796, 2.474036084054808, 0.0, 355.36545703316915, 2.3525087027248084, 1.394705056983077, 0.0, 0.029910308624236454], + linf = [1.5831869462980066, 44.20237543674303, 12.866364462538552, 0.0, 2237.8614138686, 13.066894956593798, 8.984965484247244, 0.0, 0.5226498664960756], tspan = (0.0, 0.003)) end end diff --git a/test/test_examples_2d_mhd.jl b/test/test_examples_2d_mhd.jl index 8b3e1658cd..b61d6b2df9 100644 --- a/test/test_examples_2d_mhd.jl +++ b/test/test_examples_2d_mhd.jl @@ -11,34 +11,34 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "MHD" begin @testset "elixir_mhd_alfven_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2 = [0.00011134327871126106, 5.880012991229279e-6, 5.880012991198096e-6, 8.432976937044043e-6, 1.294412985383592e-6, 1.2239229634469704e-6, 1.2239229634550153e-6, 1.8308368268137684e-6, 8.098531574321276e-7], - linf = [0.0002678911590654476, 1.6247294531021583e-5, 1.6247294531257506e-5, 2.7353154099893362e-5, 5.328011688177092e-6, 8.115802777819425e-6, 8.115802777708403e-6, 1.2104880218147263e-5, 4.179631827493849e-6]) + l2 = [0.000112085405989522, 5.908049046758828e-6, 5.908049046772597e-6, 8.493815951190774e-6, 1.297306491315482e-6, 1.2316941466595997e-6, 1.2316941466702522e-6, 1.8368487700950752e-6, 4.4903164608485606e-7], + linf = [0.0002693844074639351, 1.6289386154083596e-5, 1.6289386155082797e-5, 2.747802438275715e-5, 5.344948833307939e-6, 8.137608476288527e-6, 8.137608476843639e-6, 1.2121292898431557e-5, 2.306551665486275e-6]) end @testset "elixir_mhd_alfven_wave_mortar.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"), - l2 = [4.610806814382062e-6, 1.6898720654838787e-6, 1.6209422082463377e-6, 1.6995459258068932e-6, 1.486433605116548e-6, 1.3876687274336836e-6, 1.3412591287808953e-6, 1.715595049725351e-6, 9.813856592499105e-7], - linf = [3.5225644818837054e-5, 1.5350075873457603e-5, 1.4264708736805298e-5, 1.4421606211581506e-5, 7.744178943891455e-6, 1.0188298347757474e-5, 9.862369056090614e-6, 1.6018083168742314e-5, 5.563792310153545e-6], + l2 = [3.843336695097301e-6, 1.465357374916213e-6, 1.416088654606427e-6, 1.6806652628357514e-6, 1.4256305721007091e-6, 1.4507398918063128e-6, 1.397023085105777e-6, 1.704347968034822e-6, 8.210831113360175e-7], + linf = [3.4691159990107856e-5, 1.4347253854782305e-5, 1.39581811973849e-5, 1.4252304932890758e-5, 8.244903569987194e-6, 1.0371383968976744e-5, 1.1133207181823757e-5, 1.600372637142189e-5, 5.572707635298199e-6], tspan = (0.0, 1.0)) end @testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2 = [0.03628213711730772, 0.043011106541313036, 0.04299694609483342, 0.025746634517134097, 0.1620505609687871, 0.01745407144590488, 0.017454880786441977, 0.026880815662339404, 0.0001427849600999494], - linf = [0.2350909734343778, 0.31557861254908504, 0.30933432081715095, 0.21172448318061543, 0.9738530715930489, 0.09072863714076829, 0.09150557090076483, 0.1572617079252614, 0.004407940752443958]) + l2 = [0.03628262478842946, 0.043011555956537974, 0.042997357872762064, 0.02574553086994775, 0.1620610267679157, 0.01745438427233919, 0.017455233781000938, 0.026879482381500206, 0.00020380087569640704], + linf = [0.2351786266785303, 0.3156206778985494, 0.30941696929809714, 0.21167563519254826, 0.9688251298546073, 0.09076254289154928, 0.0916058976949814, 0.15698032974768517, 0.006131914796905391]) end @testset "elixir_mhd_orszag_tang.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), - l2 = [0.21672079513239845, 0.263644834080207, 0.31388290886723225, 0.0, 0.5122895724585859, 0.22917747391916293, 0.3430662443376359, 0.0, 0.0031626307824603], - linf = [1.2643446538250775, 0.6735720190596391, 0.8584758183037491, 0.0, 2.8094485872574957, 0.6587140271779941, 0.9599314012066199, 0.0, 0.052797557445107016], + l2 = [0.21671818699492046, 0.26364533451882755, 0.31388339884931576, 0.0, 0.5122729051056656, 0.2291828683723837, 0.34308184804821307, 0.0, 0.003210081660131638], + linf = [1.2648775120640154, 0.673836824735164, 0.8585376516993447, 0.0, 2.8148181376188246, 0.657323372510843, 0.9595541188847632, 0.0, 0.0525390951527749], tspan = (0.0, 0.09)) end @testset "elixir_mhd_orszag_tang.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), - l2 = [0.10794984293584534, 0.20177076633535806, 0.22968340939376386, 0.0, 0.29942323222232975, 0.156948176949458, 0.24306814332384577, 0.0, 0.013279172422637332], - linf = [0.5601603816370189, 0.5107264221998954, 0.6568655409846257, 0.0, 0.9903131924303785, 0.3993471648741348, 0.6835807774807094, 0.0, 0.520606845234404], + l2 = [0.10793703703273708, 0.20177434827232335, 0.2296817363329907, 0.0, 0.2994119440979309, 0.15675567943351448, 0.24281245338083307, 0.0, 0.003548624701792127], + linf = [0.560160833798496, 0.5098933680011135, 0.6566913038761877, 0.0, 0.9905839416293134, 0.39936698379939284, 0.6734754381366532, 0.0, 0.12739518123975463], tspan = (0.0, 0.06), surface_flux = flux_hll) end end From 5a7e63198a37fdcc8e468963677b731ee81b9a65 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 21:05:50 +0100 Subject: [PATCH 10/21] address many comments and simplifications --- examples/2d/elixir_mhd_alfven_wave.jl | 6 +++--- examples/2d/elixir_mhd_alfven_wave_mortar.jl | 6 +++--- examples/2d/elixir_mhd_blast_wave.jl | 6 +++--- examples/2d/elixir_mhd_ec.jl | 6 +++--- examples/2d/elixir_mhd_orszag_tang.jl | 6 +++--- examples/2d/elixir_mhd_rotor.jl | 6 +++--- src/Trixi.jl | 9 +++++---- src/callbacks_step/glm_speed.jl | 12 ++++-------- src/equations/ideal_glm_mhd_2d.jl | 3 +-- 9 files changed, 28 insertions(+), 32 deletions(-) diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl index c164608296..e3ea4ec6fe 100644 --- a/examples/2d/elixir_mhd_alfven_wave.jl +++ b/examples/2d/elixir_mhd_alfven_wave.jl @@ -44,10 +44,10 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -cfl_number = 1.0 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_alfven_wave_mortar.jl b/examples/2d/elixir_mhd_alfven_wave_mortar.jl index 4ae9a9597c..ec81016375 100644 --- a/examples/2d/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/2d/elixir_mhd_alfven_wave_mortar.jl @@ -47,10 +47,10 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -cfl_number = 1.0 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_blast_wave.jl b/examples/2d/elixir_mhd_blast_wave.jl index 2f6f58a923..8dcaaaa847 100644 --- a/examples/2d/elixir_mhd_blast_wave.jl +++ b/examples/2d/elixir_mhd_blast_wave.jl @@ -63,10 +63,10 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -cfl_number = 0.3 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 0.3 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_ec.jl b/examples/2d/elixir_mhd_ec.jl index 7a0dff6a65..d731ede9b8 100644 --- a/examples/2d/elixir_mhd_ec.jl +++ b/examples/2d/elixir_mhd_ec.jl @@ -41,10 +41,10 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -cfl_number = 1.0 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 1.0 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_orszag_tang.jl b/examples/2d/elixir_mhd_orszag_tang.jl index 4746c3c4dc..483712c1c6 100644 --- a/examples/2d/elixir_mhd_orszag_tang.jl +++ b/examples/2d/elixir_mhd_orszag_tang.jl @@ -63,10 +63,10 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -cfl_number = 0.4 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 0.4 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_rotor.jl b/examples/2d/elixir_mhd_rotor.jl index 173b76300f..b82c4c243c 100644 --- a/examples/2d/elixir_mhd_rotor.jl +++ b/examples/2d/elixir_mhd_rotor.jl @@ -64,10 +64,10 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -cfl_number = 0.35 -stepsize_callback = StepsizeCallback(cfl=cfl_number) +cfl = 0.35 +stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl_number) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/src/Trixi.jl b/src/Trixi.jl index a4fa9dc3b9..e496e6500c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -20,7 +20,8 @@ module Trixi using LinearAlgebra: dot using Printf: @printf, @sprintf, println -import DiffEqBase: ODEProblem, ODESolution, get_du, u_modified!, set_proposed_dt!, terminate! +import DiffEqBase: ODEProblem, ODESolution, get_du, u_modified!, set_proposed_dt!, terminate!, + get_proposed_dt using DiffEqCallbacks: CallbackSet, DiscreteCallback using EllipsisNotation # .. using HDF5: h5open, attrs @@ -107,9 +108,9 @@ export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integ export SemidiscretizationEulerGravity, ParametersEulerGravity, timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N! -export SummaryCallback, SteadyStateCallback, AMRCallback, StepsizeCallback, - SaveRestartCallback, SaveSolutionCallback, AnalysisCallback, AliveCallback, - TrivialCallback, GlmSpeedCallback +export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, + SaveRestartCallback, SaveSolutionCallback, AMRCallback, StepsizeCallback, + GlmSpeedCallback, TrivialCallback export load_mesh, load_time diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index ec1d2934f8..df055ecaca 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -2,10 +2,10 @@ """ GlmSpeedCallback(; glm_scale=0.5, cfl_scale=1.0) -Update the divergence cleaning wave speed c_h according to the time step -computed in StepsizeCallback for the ideal GLM-MHD equations. +Update the divergence cleaning wave speed `c_h` according to the time step +computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. """ -mutable struct GlmSpeedCallback{RealT<:Real} +struct GlmSpeedCallback{RealT<:Real} glm_scale::RealT cfl_scale::RealT end @@ -58,11 +58,7 @@ end @unpack glm_scale, cfl_scale = glm_speed_callback # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) - max_scaled_speed_for_c_h = nextfloat(zero(dt)) - for element in eachelement(solver, cache) - inv_jacobian = cache.elements.inverse_jacobian[element] - max_scaled_speed_for_c_h = max(max_scaled_speed_for_c_h, ndims(semi.equations) * inv_jacobian) - end + max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(semi.equations) c_h_deltat = cfl_scale * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) # c_h is proportional to its own time step divided by the complete MHD time step diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index d4461f7db2..02dec1163f 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -544,8 +544,7 @@ end v3 = rho_v3 / rho cf_x_direction = calc_fast_wavespeed(u, 1, equations) cf_y_direction = calc_fast_wavespeed(u, 2, equations) - cf_max = max(cf_x_direction, cf_y_direction) - + return abs(v1) + cf_x_direction, abs(v2) + cf_y_direction end From 6cfb56676a4c6e255bda54a8e25190bdf3c1f341 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 21:06:16 +0100 Subject: [PATCH 11/21] rename the new callback file --- src/callbacks_step/{glm_speed.jl => glm_speed_dg.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/callbacks_step/{glm_speed.jl => glm_speed_dg.jl} (100%) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed_dg.jl similarity index 100% rename from src/callbacks_step/glm_speed.jl rename to src/callbacks_step/glm_speed_dg.jl From 13fe408ff703219b77867ce977b1e6634fddabe0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 21:06:48 +0100 Subject: [PATCH 12/21] forgot to rename file in the include --- src/callbacks_step/callbacks_step.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index 9879fded93..f368c421d1 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -45,7 +45,7 @@ include("save_restart.jl") include("save_solution.jl") include("amr.jl") include("stepsize.jl") -include("glm_speed.jl") +include("glm_speed_dg.jl") # The `TrivialCallback` purposely does nothing: It allows to quickly disable specific callbacks From 329e216efca0952249e618165e342c0f6fa73d4e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 22:07:17 +0100 Subject: [PATCH 13/21] rename file --- src/callbacks_step/{glm_speed_dg.jl => glm_speed.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/callbacks_step/{glm_speed_dg.jl => glm_speed.jl} (100%) diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed.jl similarity index 100% rename from src/callbacks_step/glm_speed_dg.jl rename to src/callbacks_step/glm_speed.jl From 4f7fe2081e9f8e07605c1f8f0f852782c0f71ac7 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 22:08:37 +0100 Subject: [PATCH 14/21] address comments, separate computation of ch for dg into a new file, rename for consistency --- examples/2d/elixir_mhd_alfven_wave.jl | 2 +- examples/2d/elixir_mhd_alfven_wave_mortar.jl | 2 +- examples/2d/elixir_mhd_blast_wave.jl | 2 +- examples/2d/elixir_mhd_ec.jl | 2 +- examples/2d/elixir_mhd_orszag_tang.jl | 2 +- examples/2d/elixir_mhd_rotor.jl | 2 +- src/auxiliary/precompile.jl | 2 +- src/callbacks_step/callbacks_step.jl | 2 +- src/callbacks_step/glm_speed.jl | 22 +++++++++++--------- src/callbacks_step/glm_speed_dg.jl | 7 +++++++ 10 files changed, 27 insertions(+), 18 deletions(-) create mode 100644 src/callbacks_step/glm_speed_dg.jl diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl index e3ea4ec6fe..b70824049d 100644 --- a/examples/2d/elixir_mhd_alfven_wave.jl +++ b/examples/2d/elixir_mhd_alfven_wave.jl @@ -47,7 +47,7 @@ save_solution = SaveSolutionCallback(interval=10, cfl = 1.0 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_alfven_wave_mortar.jl b/examples/2d/elixir_mhd_alfven_wave_mortar.jl index ec81016375..893c50cde9 100644 --- a/examples/2d/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/2d/elixir_mhd_alfven_wave_mortar.jl @@ -50,7 +50,7 @@ save_solution = SaveSolutionCallback(interval=10, cfl = 1.0 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_blast_wave.jl b/examples/2d/elixir_mhd_blast_wave.jl index 8dcaaaa847..01428b623f 100644 --- a/examples/2d/elixir_mhd_blast_wave.jl +++ b/examples/2d/elixir_mhd_blast_wave.jl @@ -66,7 +66,7 @@ amr_callback = AMRCallback(semi, amr_controller, cfl = 0.3 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_ec.jl b/examples/2d/elixir_mhd_ec.jl index d731ede9b8..ca4fd1fea8 100644 --- a/examples/2d/elixir_mhd_ec.jl +++ b/examples/2d/elixir_mhd_ec.jl @@ -44,7 +44,7 @@ save_solution = SaveSolutionCallback(interval=10, cfl = 1.0 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_orszag_tang.jl b/examples/2d/elixir_mhd_orszag_tang.jl index 483712c1c6..16f64f021b 100644 --- a/examples/2d/elixir_mhd_orszag_tang.jl +++ b/examples/2d/elixir_mhd_orszag_tang.jl @@ -66,7 +66,7 @@ amr_callback = AMRCallback(semi, amr_controller, cfl = 0.4 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/2d/elixir_mhd_rotor.jl b/examples/2d/elixir_mhd_rotor.jl index b82c4c243c..90625b92d8 100644 --- a/examples/2d/elixir_mhd_rotor.jl +++ b/examples/2d/elixir_mhd_rotor.jl @@ -67,7 +67,7 @@ amr_callback = AMRCallback(semi, amr_controller, cfl = 0.35 stepsize_callback = StepsizeCallback(cfl=cfl) -glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl_scale=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index d5b6f7c501..ffd1f01811 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -230,7 +230,7 @@ function _precompile_manual_() @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:analysis_interval,),Tuple{Int}},Type{AliveCallback}}) for RealT in (Float64,) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:cfl,),Tuple{RealT}},Type{StepsizeCallback}}) - @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:glm_scale, :cfl_scale),Tuple{RealT,RealT}},Type{GlmSpeedCallback}}) + @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:glm_scale, :cfl),Tuple{RealT,RealT}},Type{GlmSpeedCallback}}) end @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_final_restart),Tuple{Int,Bool}},Type{SaveRestartCallback}}) @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),NamedTuple{(:interval, :save_initial_solution, :save_final_solution, :solution_variables),Tuple{Int,Bool,Bool,Symbol}},Type{SaveSolutionCallback}}) diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index f368c421d1..9879fded93 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -45,7 +45,7 @@ include("save_restart.jl") include("save_solution.jl") include("amr.jl") include("stepsize.jl") -include("glm_speed_dg.jl") +include("glm_speed.jl") # The `TrivialCallback` purposely does nothing: It allows to quickly disable specific callbacks diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index df055ecaca..78dcb3118d 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -1,20 +1,20 @@ """ - GlmSpeedCallback(; glm_scale=0.5, cfl_scale=1.0) + GlmSpeedCallback(; glm_scale=0.5, cfl=1.0) Update the divergence cleaning wave speed `c_h` according to the time step computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. """ struct GlmSpeedCallback{RealT<:Real} glm_scale::RealT - cfl_scale::RealT + cfl::RealT end function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:GlmSpeedCallback} glm_speed_callback = cb.affect! - @unpack glm_scale = glm_speed_callback - print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ")") + @unpack glm_scale, cfl = glm_speed_callback + print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, ")") end @@ -26,17 +26,18 @@ function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Af setup = [ "GLM wave speed scaling" => glm_speed_callback.glm_scale, + "Expected CFL number" => glm_speed_callback.cfl, ] summary_box(io, "GlmSpeedCallback", setup) end end -function GlmSpeedCallback(; glm_scale=0.5, cfl_scale=1.0) +function GlmSpeedCallback(; glm_scale=0.5, cfl=1.0) # when is the callback activated condition = (u, t, integrator) -> true - glm_speed_callback = GlmSpeedCallback(glm_scale, cfl_scale) + glm_speed_callback = GlmSpeedCallback(glm_scale, cfl) DiscreteCallback(condition, glm_speed_callback, save_positions=(false,false), @@ -54,15 +55,16 @@ end dt = integrator.dtcache semi = integrator.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - @unpack glm_scale, cfl_scale = glm_speed_callback + _, equations, solver, cache = mesh_equations_solver_cache(semi) + @unpack glm_scale, cfl = glm_speed_callback # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) - max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(semi.equations) - c_h_deltat = cfl_scale * 2 / (nnodes(solver) * max_scaled_speed_for_c_h) + c_h_deltat = calc_dt_for_cleaning_speed(cfl, equations, solver, cache) # c_h is proportional to its own time step divided by the complete MHD time step equations.c_h = glm_scale * c_h_deltat / dt return nothing end + +include("glm_speed_dg.jl") diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed_dg.jl new file mode 100644 index 0000000000..15f4578979 --- /dev/null +++ b/src/callbacks_step/glm_speed_dg.jl @@ -0,0 +1,7 @@ + +function calc_dt_for_cleaning_speed(cfl::Real, equations::AbstractIdealGlmMhdEquations, dg::DG, cache) +# compute time step for GLM linear advection equation with c_h=1 for the DG discretization on +# Cartesian meshes + max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(equations) + return cfl * 2 / (nnodes(dg) * max_scaled_speed_for_c_h) +end From 4d9b824725c076338876d349269bb103ef33b189 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 20 Nov 2020 22:18:55 +0100 Subject: [PATCH 15/21] forgot to update how the timestep is extracted from the integrator --- src/callbacks_step/glm_speed.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 78dcb3118d..bc7aa0e281 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -53,7 +53,7 @@ end # This method is called as callback after the StepsizeCallback during the time integration. @inline function (glm_speed_callback::GlmSpeedCallback)(integrator) - dt = integrator.dtcache + dt = get_proposed_dt(integrator) semi = integrator.p _, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack glm_scale, cfl = glm_speed_callback From 4722d531419d846011fce9aa4f93e14192df52ec Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 21 Nov 2020 00:14:09 +0100 Subject: [PATCH 16/21] updated 3d tests to use new c_h calculation --- examples/3d/elixir_mhd_alfven_wave.jl | 7 +++++-- examples/3d/elixir_mhd_alfven_wave_mortar.jl | 8 ++++++-- examples/3d/elixir_mhd_ec.jl | 10 +++++++--- examples/3d/elixir_mhd_orszag_tang.jl | 12 ++++++++---- src/callbacks_step/stepsize_dg3d.jl | 8 -------- src/equations/ideal_glm_mhd_3d.jl | 6 ------ test/test_examples_3d_mhd.jl | 20 ++++++++++---------- 7 files changed, 36 insertions(+), 35 deletions(-) diff --git a/examples/3d/elixir_mhd_alfven_wave.jl b/examples/3d/elixir_mhd_alfven_wave.jl index 987350c4bb..d77e8399f9 100644 --- a/examples/3d/elixir_mhd_alfven_wave.jl +++ b/examples/3d/elixir_mhd_alfven_wave.jl @@ -40,13 +40,16 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.4) +cfl = 1.4 +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) ############################################################################### diff --git a/examples/3d/elixir_mhd_alfven_wave_mortar.jl b/examples/3d/elixir_mhd_alfven_wave_mortar.jl index 44f87160b5..81438d7ca4 100644 --- a/examples/3d/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/3d/elixir_mhd_alfven_wave_mortar.jl @@ -43,12 +43,16 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.4) +cfl = 1.4 +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) ############################################################################### diff --git a/examples/3d/elixir_mhd_ec.jl b/examples/3d/elixir_mhd_ec.jl index a0439a0537..5403ab1335 100644 --- a/examples/3d/elixir_mhd_ec.jl +++ b/examples/3d/elixir_mhd_ec.jl @@ -41,12 +41,16 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.4) +cfl = 1.4 +stepsize_callback = StepsizeCallback(cfl=cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, + analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) ############################################################################### diff --git a/examples/3d/elixir_mhd_orszag_tang.jl b/examples/3d/elixir_mhd_orszag_tang.jl index 7296a79a2e..dd4bd78329 100644 --- a/examples/3d/elixir_mhd_orszag_tang.jl +++ b/examples/3d/elixir_mhd_orszag_tang.jl @@ -38,15 +38,19 @@ save_solution = SaveSolutionCallback(interval=10, save_final_solution=true, solution_variables=:primitive) -stepsize_callback = StepsizeCallback(cfl=1.1) +cfl = 1.1 +stepsize_callback = StepsizeCallback(cfl=cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) analysis_interval = 100 -callbacks = CallbackSet(summary_callback, - analysis_callback, +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback, + glm_speed_callback) ############################################################################### diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index 157d6d07c5..a038fabbfa 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -5,13 +5,6 @@ function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) - # FIXME: This should be implemented properly using another callback - # or something else, cf. - # https://github.com/trixi-framework/Trixi.jl/issues/257 - if equations isa AbstractIdealGlmMhdEquations - equations.c_h = zero(equations.c_h) - end - for element in eachelement(dg, cache) max_λ1 = max_λ2 = max_λ3 = zero(max_scaled_speed) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -43,4 +36,3 @@ function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, return 2 / (nnodes(dg) * max_scaled_speed) end - diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 530681b8d6..4558d9515e 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -498,12 +498,6 @@ end cf_x_direction = calc_fast_wavespeed(u, 1, equations) cf_y_direction = calc_fast_wavespeed(u, 2, equations) cf_z_direction = calc_fast_wavespeed(u, 3, equations) - cf_max = max(cf_x_direction, cf_y_direction, cf_z_direction) - - # FIXME: This should be implemented properly using another callback - # or something else, cf. - # https://github.com/trixi-framework/Trixi.jl/issues/257 - equations.c_h = max(equations.c_h, cf_max) # GLM cleaning speed = c_f return abs(v1) + cf_x_direction, abs(v2) + cf_y_direction, abs(v3) + cf_z_direction end diff --git a/test/test_examples_3d_mhd.jl b/test/test_examples_3d_mhd.jl index 6cf1f08710..870bb6f521 100644 --- a/test/test_examples_3d_mhd.jl +++ b/test/test_examples_3d_mhd.jl @@ -11,35 +11,35 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "3d") @testset "MHD" begin @testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2 = [0.01728564856700721, 0.017770558428794502, 0.01777055842879448, 0.017772303855024686, 0.07402251095395435, 0.010363317528939847, 0.010363317528939873, 0.010365251266968387, 0.00020781240461593321], - linf = [0.2648387662456203, 0.33478411844879813, 0.3347841184487984, 0.3698107321074581, 1.2338949711031062, 0.09857295382870013, 0.09857295382870068, 0.10426497383213318, 0.008020325762909617]) + l2 = [0.01728771605245583, 0.017772815631982047, 0.01777281563198207, 0.01777462671545554, 0.07415986730752369, 0.010366376257030126, 0.010366376257030144, 0.010367667726952105, 0.00025132983633341475], + linf = [0.2647949229021088, 0.33484186588194753, 0.33484186588194753, 0.3702865585024483, 1.2166553050316655, 0.10028094347112104, 0.10028094347112293, 0.10487511325516341, 0.011954627317821242]) end @testset "elixir_mhd_ec.jl with initial_condition=initial_condition_constant" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2 = [5.447723083700851e-16, 2.617532329387129e-15, 3.879181146896324e-15, 2.9617662710806352e-15, 2.684011421945918e-14, 1.371966389577962e-15, 1.3058638553193934e-15, 1.2475529686776005e-15, 1.4573570237027504e-15], - linf = [3.9968028886505635e-15, 1.4155343563970746e-14, 2.2898349882893854e-14, 1.8707257964933888e-14, 2.2737367544323206e-13, 1.2878587085651816e-14, 7.327471962526033e-15, 9.325873406851315e-15, 8.568189425129782e-15], + l2 = [4.270231310667203e-16, 2.4381208042014784e-15, 5.345107673575357e-15, 3.00313882171883e-15, 1.7772703118758417e-14, 1.0340110783830874e-15, 1.1779095371939702e-15, 9.961878521814573e-16, 8.1201730630719145e-16], + linf = [2.4424906541753444e-15, 2.881028748902281e-14, 2.4646951146678475e-14, 2.3092638912203256e-14, 2.3447910280083306e-13, 1.7763568394002505e-14, 1.0436096431476471e-14, 2.042810365310288e-14, 7.057203733035201e-15], atol = 1000*eps(), initial_condition=initial_condition_constant) end @testset "elixir_mhd_alfven_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2 = [0.0038734936431132137, 0.0090374494786511, 0.004173235737356234, 0.011605032955462995, 0.006247939160992442, 0.009226153647367723, 0.003460561919838917, 0.011684984122517875, 0.002201096128779557], - linf = [0.012630239144833078, 0.03265663914470475, 0.01291368559476544, 0.04444329719972474, 0.027974787995271644, 0.03453507441133391, 0.01022512252706076, 0.04498328542685666, 0.009861640501698054]) + l2 = [0.003872168982949962, 0.009234032031394519, 0.004166855440339939, 0.01160194799648245, 0.006989139550604536, 0.009258319391439123, 0.004972194549900642, 0.011701012705405973, 0.003796916988492686], + linf = [0.012897588621946787, 0.03317614952204696, 0.011594778657069232, 0.04437495747712844, 0.025695368105128846, 0.03336443221122043, 0.01000776455308272, 0.044820519469304224, 0.01606741448289395]) end @testset "elixir_mhd_alfven_wave_mortar.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"), - l2 = [0.0021484230386442484, 0.006826974752320685, 0.0030655172071723945, 0.008736026827927016, 0.005159744590177398, 0.007158578094946651, 0.0028289753083129443, 0.008815093448251163, 0.0022268140033108634], - linf = [0.013181862339541217, 0.05433825341881836, 0.02085489712458487, 0.05947341475014939, 0.03171319598340849, 0.05439174028778537, 0.017933784069765202, 0.06036923226300564, 0.012816582427655984], + l2 = [0.0021525855794505733, 0.006723734997085876, 0.0030402457081513405, 0.008711037977296627, 0.005507951439477701, 0.007204815741889562, 0.0029726037822043116, 0.008823488862283192, 0.003551010115322961], + linf = [0.013133814395385413, 0.05541886629871225, 0.020708270215062375, 0.05932492264609865, 0.03448645356854563, 0.05319127664438378, 0.01770393941703363, 0.060904784306590354, 0.020639279160633826], tspan = (0.0, 0.25)) end @testset "elixir_mhd_orszag_tang.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), - l2 = [0.0043911605686459704, 0.041447356655813394, 0.04150129977673066, 0.041503536105048366, 0.03693119824334232, 0.02112559892198415, 0.03295606821170978, 0.03296235617354949, 6.360380099124839e-6], - linf = [0.01789383890583951, 0.0848675349327856, 0.08910602912506882, 0.08491879187575965, 0.10444596146695251, 0.05381953967385888, 0.08847783436169666, 0.07784630781912419, 8.236241065117021e-5], + l2 = [0.004391160574957181, 0.041447356654791025, 0.0415012997738177, 0.04150353610875128, 0.036931198261162486, 0.021125598921156254, 0.03295606813944825, 0.032962356174273615, 6.036824268879037e-6], + linf = [0.017893840542084788, 0.08486753617328145, 0.08910603135148028, 0.08491878976510482, 0.10444599411185623, 0.05381959529402853, 0.08847807738420276, 0.07784629959770338, 7.694931670095901e-5], tspan = (0.0, 0.06)) end end From 145afafb9d328f1a90ca5ed01174b1a662256f56 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 23 Nov 2020 07:42:10 +0100 Subject: [PATCH 17/21] slight changes to make glm callback structure match the other callback structures! --- src/callbacks_step/glm_speed.jl | 4 ++-- src/callbacks_step/glm_speed_dg.jl | 5 ++++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index bc7aa0e281..2f402e33de 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -55,11 +55,11 @@ end dt = get_proposed_dt(integrator) semi = integrator.p - _, equations, solver, cache = mesh_equations_solver_cache(semi) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack glm_scale, cfl = glm_speed_callback # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR) - c_h_deltat = calc_dt_for_cleaning_speed(cfl, equations, solver, cache) + c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache) # c_h is proportional to its own time step divided by the complete MHD time step equations.c_h = glm_scale * c_h_deltat / dt diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed_dg.jl index 15f4578979..207c1a1aae 100644 --- a/src/callbacks_step/glm_speed_dg.jl +++ b/src/callbacks_step/glm_speed_dg.jl @@ -1,7 +1,10 @@ -function calc_dt_for_cleaning_speed(cfl::Real, equations::AbstractIdealGlmMhdEquations, dg::DG, cache) +function calc_dt_for_cleaning_speed(cfl::Real, mesh, + equations::AbstractIdealGlmMhdEquations, dg::DG, cache) # compute time step for GLM linear advection equation with c_h=1 for the DG discretization on # Cartesian meshes max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(equations) + # OBS! This depends on the implementation details of the StepsizeCallback and needs to be adapted + # as well if that callback changes. return cfl * 2 / (nnodes(dg) * max_scaled_speed_for_c_h) end From 1edaabd50d66f3f1450f3dbe8d9b86502e0ef9d5 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 23 Nov 2020 08:00:05 +0100 Subject: [PATCH 18/21] added assert error and more documentation for the GLM speed callback --- src/callbacks_step/glm_speed.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 2f402e33de..280c759dfa 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -4,6 +4,10 @@ Update the divergence cleaning wave speed `c_h` according to the time step computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. +The `cfl` number should be set to the same value as for the time step size calculation. The +`glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD +solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale` = 0 +deactivates the divergence cleaning. """ struct GlmSpeedCallback{RealT<:Real} glm_scale::RealT @@ -39,6 +43,8 @@ function GlmSpeedCallback(; glm_scale=0.5, cfl=1.0) glm_speed_callback = GlmSpeedCallback(glm_scale, cfl) + @assert(0<=glm_speed_callback.glm_scale<=1,"glm_scale must be between 0 and 1") + DiscreteCallback(condition, glm_speed_callback, save_positions=(false,false), initialize=initialize!) From 211cd60e7ee2a47f42f96d4c9aebf4fda337da8d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 23 Nov 2020 08:38:34 +0100 Subject: [PATCH 19/21] added a check to ensure that the c_h is set to something in the ideal MHD constructor --- src/equations/ideal_glm_mhd_2d.jl | 7 ++++--- src/equations/ideal_glm_mhd_3d.jl | 5 +++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 02dec1163f..00f3c81db1 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -9,8 +9,9 @@ mutable struct IdealGlmMhdEquations2D{RealT<:Real} <: AbstractIdealGlmMhdEquatio c_h::RealT # GLM cleaning speed end -function IdealGlmMhdEquations2D(gamma) - IdealGlmMhdEquations2D(gamma, zero(gamma)) +function IdealGlmMhdEquations2D(gamma; initial_c_h=convert(typeof(gamma), NaN)) + # Use `promote` to ensure that `gamma` and `initial_c_h` have the same type + IdealGlmMhdEquations2D(promote(gamma, initial_c_h)...) end @@ -544,7 +545,7 @@ end v3 = rho_v3 / rho cf_x_direction = calc_fast_wavespeed(u, 1, equations) cf_y_direction = calc_fast_wavespeed(u, 2, equations) - + return abs(v1) + cf_x_direction, abs(v2) + cf_y_direction end diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 4558d9515e..09087e5b00 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -9,8 +9,9 @@ mutable struct IdealGlmMhdEquations3D{RealT<:Real} <: AbstractIdealGlmMhdEquatio c_h::RealT # GLM cleaning speed end -function IdealGlmMhdEquations3D(gamma) - IdealGlmMhdEquations3D(gamma, zero(gamma)) +function IdealGlmMhdEquations3D(gamma; initial_c_h=convert(typeof(gamma), NaN)) + # Use `promote` to ensure that `gamma` and `initial_c_h` have the same type + IdealGlmMhdEquations3D(promote(gamma, initial_c_h)...) end From e84ecd0ba8755aa12e5c9639147e312fd6878a58 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 23 Nov 2020 08:39:13 +0100 Subject: [PATCH 20/21] moved assert error to be before new memory is used --- src/callbacks_step/glm_speed.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 280c759dfa..2d119f1e47 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -37,13 +37,13 @@ function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Af end -function GlmSpeedCallback(; glm_scale=0.5, cfl=1.0) +function GlmSpeedCallback(; glm_scale=0.5, cfl) # when is the callback activated condition = (u, t, integrator) -> true - glm_speed_callback = GlmSpeedCallback(glm_scale, cfl) + @assert 0 <= glm_speed_callback.glm_scale <= 1 "glm_scale must be between 0 and 1" - @assert(0<=glm_speed_callback.glm_scale<=1,"glm_scale must be between 0 and 1") + glm_speed_callback = GlmSpeedCallback(glm_scale, cfl) DiscreteCallback(condition, glm_speed_callback, save_positions=(false,false), From 01b8361d099f9fc50e7b0f3f155f67310c207e7e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 23 Nov 2020 09:05:55 +0100 Subject: [PATCH 21/21] fixed error in assert and updated docstring --- src/callbacks_step/glm_speed.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 2d119f1e47..c270a0454f 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -1,6 +1,6 @@ """ - GlmSpeedCallback(; glm_scale=0.5, cfl=1.0) + GlmSpeedCallback(; glm_scale=0.5, cfl) Update the divergence cleaning wave speed `c_h` according to the time step computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. @@ -41,7 +41,7 @@ function GlmSpeedCallback(; glm_scale=0.5, cfl) # when is the callback activated condition = (u, t, integrator) -> true - @assert 0 <= glm_speed_callback.glm_scale <= 1 "glm_scale must be between 0 and 1" + @assert 0 <= glm_scale <= 1 "glm_scale must be between 0 and 1" glm_speed_callback = GlmSpeedCallback(glm_scale, cfl)