NaN check callback

examples/structured_2d_dgsem/elixir_hypdiff_nonperiodic.jl

@@ -64,11 +64,14 @@ save_solution = SaveSolutionCallback(interval = 4000,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 stepsize_callback = StepsizeCallback(cfl = 1.9)
 
 callbacks = CallbackSet(summary_callback, steady_state_callback,
                         analysis_callback, alive_callback,
                         save_solution,
+                        nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl

@@ -70,9 +70,11 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 1.75)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback, steady_state_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_advection_timeintegration.jl

@@ -49,9 +49,11 @@ amr_callback = AMRCallback(semi, amr_controller,
 
 stepsize_callback = StepsizeCallback(cfl = 1.6)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         amr_callback, stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl

@@ -77,9 +77,11 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 0.3)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl

@@ -75,9 +75,11 @@ save_restart = SaveRestartCallback(interval = 1000,
 
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        stepsize_callback,
+                        stepsize_callback, nan_callback,
                         save_restart, save_solution)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -76,9 +76,11 @@ save_solution = SaveSolutionCallback(interval = 1000,
 
 stepsize_callback = StepsizeCallback(cfl = 0.6)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        stepsize_callback,
+                        stepsize_callback, nan_callback,
                         save_solution)
 ###############################################################################
 # run the simulation

examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl

@@ -74,9 +74,11 @@ save_solution = SaveSolutionCallback(dt = 0.1,
 
 stepsize_callback = StepsizeCallback(cfl = 0.6)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl

@@ -123,11 +123,14 @@ save_solution = SaveSolutionCallback(interval = 600,
 
 stepsize_callback = StepsizeCallback(cfl = 0.5)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
                         alive_callback,
                         save_solution,
-                        stepsize_callback)
+                        stepsize_callback,
+                        nan_callback)
 
 ###############################################################################
 # run the simulation

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl

@@ -124,11 +124,14 @@ save_solution = SaveSolutionCallback(interval = 300,
 
 stepsize_callback = StepsizeCallback(cfl = 0.9)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
                         alive_callback,
                         save_solution,
-                        stepsize_callback)
+                        stepsize_callback,
+                        nan_callback)
 
 ###############################################################################
 # run the simulation

examples/tree_2d_dgsem/elixir_hypdiff_lax_friedrichs.jl

@@ -75,9 +75,11 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 1.2)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback, steady_state_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl

@@ -88,13 +88,16 @@ save_solution = SaveSolutionCallback(interval = 100,
 cfl = 0.4
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
                         alive_callback,
                         save_solution,
                         stepsize_callback,
+                        nan_callback,
                         glm_speed_callback)
 
 ###############################################################################

examples/tree_3d_dgsem/elixir_hypdiff_lax_friedrichs.jl

@@ -79,9 +79,11 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 2.4)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback, steady_state_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_3d_dgsem/elixir_hypdiff_nonperiodic.jl

@@ -52,9 +52,11 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 1.8)
 
+nan_callback = NaNCallback(analysis_interval = analysis_interval)
+
 callbacks = CallbackSet(summary_callback, steady_state_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        save_solution, nan_callback,
                         stepsize_callback)
 
 ###############################################################################

src/Trixi.jl

@@ -259,7 +259,7 @@ export SemidiscretizationCoupled
 export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback,
        AveragingCallback,
-       AMRCallback, StepsizeCallback,
+       AMRCallback, StepsizeCallback, NaNCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled
 

src/callbacks_step/callbacks_step.jl

@@ -53,6 +53,7 @@ include("summary.jl")
 include("steady_state.jl")
 include("analysis.jl")
 include("alive.jl")
+include("nan.jl")
 include("save_restart.jl")
 include("save_solution.jl")
 include("time_series.jl")

src/callbacks_step/nan.jl

@@ -0,0 +1,85 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+"""
+    NaNCallback(analysis_interval=0, nan_interval=analysis_interval÷10)
+
+Callback checking for NaNs in the solution vector `u` every `nan_interval`.
+This should be used when using `Trixi.solve` to solve the ODEs as the custom 
+integrators do not come with the `unstable_check` of the DifferentialEquations.jl integrators.
+If `analysis_interval ≂̸ 0`, the output is omitted every `analysis_interval` time steps.
+"""
+mutable struct NaNCallback
+    nan_interval::Int
+    analysis_interval::Int
+end
+
+function NaNCallback(; analysis_interval = 0,
+                     nan_interval = analysis_interval ÷ 10)
+    nan_callback = NaNCallback(nan_interval, analysis_interval)
+
+    DiscreteCallback(nan_callback, nan_callback, # the first one is the condition, the second the affect!
+                     save_positions = (false, false),
+                     initialize = initialize!)
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:NaNCallback})
+    @nospecialize cb # reduce precompilation time
+
+    nan_callback = cb.affect!
+    print(io, "NaNCallback(nan_interval=", nan_callback.nan_interval, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, <:NaNCallback})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        nan_callback = cb.affect!
+
+        setup = [
+            "interval" => nan_callback.nan_interval,
+        ]
+        summary_box(io, "NaNCallback", setup)
+    end
+end
+
+function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t,
+                     integrator) where {Condition, Affect! <: NaNCallback}
+    nan_callback = cb.affect!
+    return nothing
+end
+
+# this method is called to determine whether the callback should be activated
+function (nan_callback::NaNCallback)(u, t, integrator)
+    @unpack nan_interval, analysis_interval = nan_callback
+
+    # With error-based step size control, some steps can be rejected. Thus,
+    #   `integrator.iter >= integrator.stats.naccept`
+    #    (total #steps)       (#accepted steps)
+    # We need to check the number of accepted steps since callbacks are not
+    # activated after a rejected step.
+    return nan_interval > 0 && ((integrator.stats.naccept % nan_interval == 0 &&
+             !(integrator.stats.naccept == 0 && integrator.iter > 0) &&
+             (analysis_interval == 0 ||
+              integrator.stats.naccept % analysis_interval != 0)) ||
+            isfinished(integrator))
+end
+
+# this method is called when the callback is activated
+function (nan_callback::NaNCallback)(integrator)
+    @trixi_timeit timer() "NaNCallback" begin
+        if any(isnan, integrator.u)
+            error("NaN detected in the solution vector `u`` at time $(integrator.t),
+                  timestep $(integrator.iter)")
+        end
+    end
+    return nothing
+end
+end # @muladd

src/time_integration/methods_2N.jl

@@ -80,7 +80,7 @@ end
 # which are used in Trixi.jl.
 mutable struct SimpleIntegrator2N{RealT <: Real, uType, Params, Sol, F, Alg,
                                   SimpleIntegrator2NOptions}
-    u::uType #
+    u::uType
     du::uType
     u_tmp::uType
     t::RealT

src/time_integration/methods_3Sstar.jl

@@ -146,7 +146,7 @@ end
 
 mutable struct SimpleIntegrator3Sstar{RealT <: Real, uType, Params, Sol, F, Alg,
                                       SimpleIntegrator3SstarOptions}
-    u::uType #
+    u::uType
     du::uType
     u_tmp1::uType
     u_tmp2::uType