Outsource saving of bounds check errors

src/callbacks_stage/subcell_bounds_check.jl

@@ -42,23 +42,35 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
                   (stage == length(alg.c)) &&
                   (iter % callback.interval == 0 || integrator.finalstep)
     @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
-                                                      solver.volume_integral, t,
-                                                      iter + 1,
-                                                      callback.output_directory,
-                                                      save_errors)
+                                                      solver.volume_integral)
+
+    if save_errors
+        @trixi_timeit timer() "save_errors" save_bounds_check_errors(callback.output_directory,
+                                                                     u, t, iter + 1,
+                                                                     equations,
+                                                                     solver.volume_integral)
+    end
+end
+
+@inline function check_bounds(u, mesh, equations, solver, cache,
+                              volume_integral::AbstractVolumeIntegral)
+    return nothing
+end
+
+@inline function check_bounds(u, mesh, equations, solver, cache,
+                              volume_integral::VolumeIntegralSubcellLimiting)
+    check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter)
 end
 
-function check_bounds(u, mesh, equations, solver, cache,
-                      volume_integral::AbstractVolumeIntegral, t, iter,
-                      output_directory, save_errors)
+@inline function save_bounds_check_errors(output_directory, u, t, iter, equations,
+                                          volume_integral::AbstractVolumeIntegral)
     return nothing
 end
 
-function check_bounds(u, mesh, equations, solver, cache,
-                      volume_integral::VolumeIntegralSubcellLimiting, t, iter,
-                      output_directory, save_errors)
-    check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter,
-                 output_directory, save_errors)
+@inline function save_bounds_check_errors(output_directory, u, t, iter, equations,
+                                          volume_integral::VolumeIntegralSubcellLimiting)
+    save_bounds_check_errors(output_directory, u, t, iter, equations,
+                             volume_integral.limiter)
 end
 
 function init_callback(callback::BoundsCheckCallback, semi)

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -6,8 +6,7 @@
 #! format: noindent
 
 @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
-                              limiter::SubcellLimiterIDP,
-                              time, iter, output_directory, save_errors)
+                              limiter::SubcellLimiterIDP)
     (; local_minmax, positivity) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache
@@ -85,40 +84,48 @@
                                            idp_bounds_delta_local[key][stride_size])
     end
 
-    if save_errors
-        # Print to output file
-        open("$output_directory/deviations.txt", "a") do f
-            print(f, iter, ", ", time)
-            if local_minmax
-                for v in limiter.local_minmax_variables_cons
-                    v_string = string(v)
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(v_string, "_min")][stride_size],
-                          ", ",
-                          idp_bounds_delta_local[Symbol(v_string, "_max")][stride_size])
-                end
+    return nothing
+end
+
+@inline function save_bounds_check_errors(output_directory, u, time, iter, equations,
+                                          limiter::SubcellLimiterIDP)
+    (; local_minmax, positivity) = limiter
+    (; idp_bounds_delta_local) = limiter.cache
+
+    stride_size = div(128, sizeof(eltype(u))) # = n
+
+    # Print errors to output file
+    open("$output_directory/deviations.txt", "a") do f
+        print(f, iter, ", ", time)
+        if local_minmax
+            for v in limiter.local_minmax_variables_cons
+                v_string = string(v)
+                print(f, ", ",
+                      idp_bounds_delta_local[Symbol(v_string, "_min")][stride_size],
+                      ", ",
+                      idp_bounds_delta_local[Symbol(v_string, "_max")][stride_size])
             end
-            if positivity
-                for v in limiter.positivity_variables_cons
-                    if v in limiter.local_minmax_variables_cons
-                        continue
-                    end
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size])
-                end
-                for variable in limiter.positivity_variables_nonlinear
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size])
+        end
+        if positivity
+            for v in limiter.positivity_variables_cons
+                if v in limiter.local_minmax_variables_cons
+                    continue
                 end
+                print(f, ", ",
+                      idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size])
             end
-            println(f)
-        end
-        # Reset local maximum deviations
-        for (key, _) in idp_bounds_delta_local
-            for i in 1:Threads.nthreads()
-                idp_bounds_delta_local[key][stride_size * i] = zero(eltype(idp_bounds_delta_local[key][stride_size]))
+            for variable in limiter.positivity_variables_nonlinear
+                print(f, ", ",
+                      idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size])
             end
         end
+        println(f)
+    end
+    # Reset local maximum deviations
+    for (key, _) in idp_bounds_delta_local
+        for i in 1:Threads.nthreads()
+            idp_bounds_delta_local[key][stride_size * i] = zero(eltype(idp_bounds_delta_local[key][stride_size]))
+        end
     end
 
     return nothing