Stabilize simulations to avoid differences between CI and local runs

docs/literate/src/files/subcell_shock_capturing.jl

@@ -51,7 +51,11 @@
 # The Newton-bisection algorithm is an iterative method and requires some parameters.
 # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and
 # relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are
-# sufficient in most cases and therefore used as default. Additionally, there is the parameter
+# sufficient in most cases and therefore used as default. If the implemented bounds checking
+# functionality indicates problems with the limiting (see [below](@ref subcell_bounds_check))
+# the Newton method with the chosen parameters might not manage to converge. If so, adapting
+# the mentioned parameters helps fix that.
+# Additionally, there is the parameter
 # `gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation
 # of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`,
 # as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.]
@@ -244,7 +248,7 @@ plot(sol)
 # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116)
 
 
-# ## Bounds checking
+# ## [Bounds checking](@id subcell_bounds_check)
 # Subcell limiting is based on the fulfillment of target bounds - either global or local.
 # Although the implementation works and has been thoroughly tested, there are some cases where
 # these bounds are not met.

examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -51,7 +51,9 @@ limiter_idp = SubcellLimiterIDP(equations, basis;
                                 local_twosided_variables_cons = ["rho"],
                                 local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
                                                                        min)],
-                                max_iterations_newton = 40, # Default value of 10 iterations is too low to fulfill bounds.
+                                # Default parameters are not sufficient to fulfill bounds properly.
+                                max_iterations_newton = 40,
+                                newton_tolerances = (1.0e-13, 1.0e-15),
                                 positivity_variables_cons = [],
                                 positivity_variables_nonlinear = [])
 # Variables for global limiting (`positivity_variables_cons` and
@@ -96,7 +98,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.7)
+stepsize_callback = StepsizeCallback(cfl = 0.6)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl

@@ -43,7 +43,10 @@ basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 local_twosided_variables_cons = ["rho"],
                                 local_onesided_variables_nonlinear = [(Trixi.entropy_math,
-                                                                       max)])
+                                                                       max)],
+                                # Default parameters are not sufficient to fulfill bounds properly.
+                                max_iterations_newton = 70,
+                                newton_tolerances = (1.0e-13, 1.0e-14))
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -44,7 +44,11 @@ basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 local_twosided_variables_cons = ["rho"],
                                 local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
-                                                                       min)])
+                                                                       min)],
+                                positivity_variables_nonlinear = [pressure],
+                                # Default parameters are not sufficient to fulfill bounds properly.
+                                max_iterations_newton = 60,
+                                newton_tolerances = (1.0e-13, 1.0e-15))
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
@@ -53,7 +57,7 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (-2.0, -2.0)
 coordinates_max = (2.0, 2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 3,
+                initial_refinement_level = 5,
                 n_cells_max = 100_000)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
@@ -76,7 +80,7 @@ save_solution = SaveSolutionCallback(interval = 1000,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.6)
+stepsize_callback = StepsizeCallback(cfl = 0.4)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl

@@ -121,7 +121,7 @@ save_solution = SaveSolutionCallback(interval = 600,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.5)
+stepsize_callback = StepsizeCallback(cfl = 0.4)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -329,9 +329,6 @@ end
             Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) +
                  min(0, val_flux2_local) + min(0, val_flux2_local_jp1)
 
-            Qp = max(0, (var_max[i, j, element] - var) / dt)
-            Qm = min(0, (var_min[i, j, element] - var) / dt)
-
             Pp = inverse_jacobian * Pp
             Pm = inverse_jacobian * Pm
 

test/test_structured_2d.jl

@@ -647,16 +647,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
                         l2=[
-                            0.6337774834710513,
-                            0.30377119245852724,
-                            0.3111372568571772,
-                            1.2976221893997268,
+                            0.6403528328480915,
+                            0.3068073114438902,
+                            0.3140151910019577,
+                            1.2977732581465693,
                         ],
                         linf=[
-                            2.2064877103138207,
-                            1.541067099687334,
-                            1.5487587769900337,
-                            6.271271639873466,
+                            2.239791987419344,
+                            1.5580885989144924,
+                            1.5392923786831547,
+                            6.2729281824590855,
                         ],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations
@@ -677,16 +677,16 @@ end
                         local_twosided_variables_cons=[],
                         local_onesided_variables_nonlinear=[],
                         l2=[
-                            0.7869912572385168,
-                            0.39170886758882073,
-                            0.39613257454431977,
-                            1.2951760266455101,
+                            0.7979084213982606,
+                            0.3980284851419719,
+                            0.4021949448633982,
+                            1.2956482394747346,
                         ],
                         linf=[
-                            5.156044534854053,
-                            3.6261667239538986,
-                            3.1807681416546085,
-                            6.3028422220287235,
+                            5.477809925838038,
+                            3.7793130706228273,
+                            3.2838862964081637,
+                            6.316943647948965,
                         ],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations

test/test_tree_2d_euler.jl

@@ -349,16 +349,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"),
                         l2=[
-                            0.3221177942225801,
-                            0.1798478357478982,
-                            0.1798364616438908,
-                            0.6136884131056267,
+                            0.3221078812528291,
+                            0.17985175694043076,
+                            0.17983453493705628,
+                            0.6136916718599121,
                         ],
                         linf=[
-                            1.343766644801395,
-                            1.1749593109683463,
-                            1.1747613085307178,
-                            2.4216006041018785,
+                            1.343237509126809,
+                            1.1747101056222315,
+                            1.174585608472406,
+                            2.4216027326405487,
                         ],
                         tspan=(0.0, 0.5),
                         initial_refinement_level=4,
@@ -404,32 +404,33 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
                         l2=[
-                            0.41444427153173785,
-                            0.1460669409661223,
-                            0.14606693069201596,
-                            0.6168046457461059,
+                            0.4227549115123529,
+                            0.14825759222652649,
+                            0.14825759222682933,
+                            0.6164668313131949,
                         ],
                         linf=[
-                            1.5720584643579567,
-                            0.7946656826861964,
-                            0.7946656525739751,
-                            6.455520291414711,
+                            1.6391908143728386,
+                            0.8344433355906021,
+                            0.8344433355966195,
+                            6.450305752671201,
                         ],
                         tspan=(0.0, 1.0),
                         initial_refinement_level=4,
                         coverage_override=(maxiters = 6,),
                         save_errors=true)
     lines = readlines(joinpath("out", "deviations.txt"))
-    @test lines[1] == "# iter, simu_time, rho_min, rho_max, entropy_guermond_etal_min"
+    @test lines[1] ==
+          "# iter, simu_time, rho_min, rho_max, entropy_guermond_etal_min, pressure_min"
     cmd = string(Base.julia_cmd())
     coverage = occursin("--code-coverage", cmd) &&
                !occursin("--code-coverage=none", cmd)
     if coverage
         # Run with coverage takes 6 time steps.
         @test startswith(lines[end], "6")
     else
-        # Run without coverage takes 89 time steps.
-        @test startswith(lines[end], "89")
+        # Run without coverage takes 138 time steps.
+        @test startswith(lines[end], "138")
     end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

test/test_tree_2d_eulermulti.jl

@@ -99,18 +99,18 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"),
                             l2=[
-                                73.10860950390489,
-                                1.4599090197303102,
-                                57176.23978426408,
-                                0.17812910616624406,
-                                0.010123079422717837,
+                                73.41054363926742,
+                                1.5072038797716156,
+                                57405.58964098063,
+                                0.17877099207437225,
+                                0.010085388785440972,
                             ],
                             linf=[
-                                214.50568817511956,
-                                25.40392579616452,
-                                152862.41011222568,
-                                0.564195553101797,
-                                0.0956331651771212,
+                                213.59140793740318,
+                                24.57625853486584,
+                                152498.21319871658,
+                                0.5911106543157919,
+                                0.09936092838440383,
                             ],
                             initial_refinement_level=3,
                             tspan=(0.0, 0.001))

test/test_tree_2d_mhd.jl

@@ -355,8 +355,8 @@ end
                             6.1013422157115546e-03,
                         ],
                         tspan=(0.0, 0.003))
-    # Ensure that we do not have excessive memory allocations 
-    # (e.g., from type instabilities) 
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
     let
         t = sol.t[end]
         u_ode = sol.u[end]