Merge pull request trixi-framework#353 from trixi-framework/fix_calc_ch

Fix calc ch
ranocha · Nov 23, 2020 · 79c1a64 · 79c1a64
2 parents eaa53dd + 01b8361
commit 79c1a64
Show file tree

Hide file tree

Showing 22 changed files with 205 additions and 102 deletions.
diff --git a/examples/2d/elixir_mhd_alfven_wave.jl b/examples/2d/elixir_mhd_alfven_wave.jl
@@ -36,22 +36,25 @@ 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,
                                      save_initial_solution=true,
                                      save_final_solution=true,
                                      solution_variables=:primitive)
 
-stepsize_callback = StepsizeCallback(cfl=1.0)
+cfl = 1.0
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
 
-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/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 = 1.0
+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)
 
 ###############################################################################
 # run the simulation

diff --git 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 = 0.3
+stepsize_callback = StepsizeCallback(cfl=cfl)
 
-callbacks = CallbackSet(summary_callback, 
-                        analysis_callback, 
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+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
@@ -41,13 +41,17 @@ save_solution = SaveSolutionCallback(interval=10,
                                      save_final_solution=true,
                                      solution_variables=:primitive)
 
-stepsize_callback = StepsizeCallback(cfl=1.0)
+cfl = 1.0
+stepsize_callback = StepsizeCallback(cfl=cfl)
 
-callbacks = CallbackSet(summary_callback, 
-                        analysis_callback, 
+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/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 = 0.4
+stepsize_callback = StepsizeCallback(cfl=cfl)
 
-callbacks = CallbackSet(summary_callback, 
-                        analysis_callback, 
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+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
@@ -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 = 0.35
+stepsize_callback = StepsizeCallback(cfl=cfl)
 
-callbacks = CallbackSet(summary_callback, 
-                        analysis_callback, 
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+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/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
@@ -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
@@ -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
@@ -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/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, attributes
@@ -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
+export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
+       SaveRestartCallback, SaveSolutionCallback, AMRCallback, StepsizeCallback,
+       GlmSpeedCallback, TrivialCallback
 
 export load_mesh, load_time
 

diff --git 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, :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}})
@@ -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
@@ -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
@@ -0,0 +1,76 @@
+
+"""
+    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.
+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
+  cfl::RealT
+end
+
+
+function Base.show(io::IO, cb::DiscreteCallback{Condition,Affect!}) where {Condition, Affect!<:GlmSpeedCallback}
+  glm_speed_callback = cb.affect!
+  @unpack glm_scale, cfl = glm_speed_callback
+  print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, ")")
+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,
+             "Expected CFL number" => glm_speed_callback.cfl,
+            ]
+    summary_box(io, "GlmSpeedCallback", setup)
+  end
+end
+
+
+function GlmSpeedCallback(; glm_scale=0.5, cfl)
+  # when is the callback activated
+  condition = (u, t, integrator) -> true
+
+  @assert 0 <= 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),
+                   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 = get_proposed_dt(integrator)
+  semi = integrator.p
+  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, 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
+
+  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
@@ -0,0 +1,10 @@
+
+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
diff --git 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)