diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 8b5399c23a..f47cdb6996 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/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. diff --git a/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 5c11a7d15a..ce603359be 100644 --- a/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/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, diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl index 00d3c69f2e..21a89fdcbe 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl +++ b/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) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 2089d35397..1dc0586ec7 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/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, diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl index b2d49ecbd4..9bf310d753 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl +++ b/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, diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 5dae5798a8..df24a58885 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/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 diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 9814b73dd4..ba2565e0dd 100644 --- a/test/test_structured_2d.jl +++ b/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 @@ -665,8 +665,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 10000 @@ -681,16 +681,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 @@ -699,8 +699,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 10000 @@ -876,8 +876,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 85ae56152d..6685c289a9 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -231,8 +231,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -267,8 +267,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -357,16 +357,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, @@ -377,8 +377,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -416,23 +416,24 @@ 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) @@ -440,8 +441,8 @@ end # 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) @@ -449,8 +450,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -672,8 +673,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index ad6ae04f83..d6dc3ef94f 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -91,8 +91,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -103,18 +103,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)) @@ -124,8 +124,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index 79c79dfb94..ed6b7f6b52 100644 --- a/test/test_tree_2d_mhd.jl +++ b/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]