From d772bf88d3c3693c3228e68b65eb1df31ea30a98 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Mon, 13 Nov 2023 06:08:09 +0100 Subject: [PATCH 1/8] CompatHelper: bump compat for Trixi to 0.6 for package benchmark, (keep existing compat) (#1735) Co-authored-by: CompatHelper Julia --- benchmark/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/Project.toml b/benchmark/Project.toml index fb985572532..e94144cfd15 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -8,4 +8,4 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" BenchmarkTools = "0.5, 0.7, 1.0" OrdinaryDiffEq = "5.65, 6" PkgBenchmark = "0.2.10" -Trixi = "0.4, 0.5" +Trixi = "0.4, 0.5, 0.6" From 7bb3b463a660edb5754fcf917ff8812f6791716a Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Mon, 13 Nov 2023 14:57:06 +0100 Subject: [PATCH 2/8] Add two-sided Zalesak-type IDP subcell limiting (#1648) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add two-sided limiting for conservative variables * Fix visualization routines * Add bounds calculation for BoundaryConditionDirichlet * Reduce cfl in elixir * Fix test * Add comment about subcell limiting with non-conforming meshes * Remove subcell visualization * Fix last commit * Remove @unpack * Add comment to `News.md` * Fix source for sedov blast setup; Formatting * Reduce allocations * Replace construction of Symbols * Add bounds check for local limiting * Implement suggestions * Fix format * Add subcell allocation tests; Add changes to minmax limiter * Undo changes in elixirs * Implement suggestions * Skip positivity limiting if local limiting is more restrictive * Reduce allocations * Pass variables as strings instead of ints * Add `_nonperiodic` to elixir name * Fix unit test * Implement suggestions * Add missing comma in export of bounds check deviation * Implement suggestions --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Andrés Rueda-Ramírez --- NEWS.md | 3 +- ...euler_blast_wave_sc_subcell_nonperiodic.jl | 93 +++++++++ ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 91 +++++++++ .../elixir_euler_shockcapturing_subcell.jl | 2 +- ...ck_bubble_shockcapturing_subcell_minmax.jl | 142 ++++++++++++++ ...ubble_shockcapturing_subcell_positivity.jl | 5 +- .../elixir_mhd_shockcapturing_subcell.jl | 2 +- src/callbacks_stage/subcell_bounds_check.jl | 24 ++- .../subcell_bounds_check_2d.jl | 38 +++- src/equations/equations.jl | 20 +- src/solvers/dg.jl | 6 + .../dgsem_tree/dg_2d_subcell_limiters.jl | 17 ++ src/solvers/dgsem_tree/subcell_limiters.jl | 66 +++++-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 178 ++++++++++++++++++ test/test_tree_2d_euler.jl | 56 ++++++ test/test_tree_2d_eulermulti.jl | 29 +++ test/test_unit.jl | 2 +- 17 files changed, 744 insertions(+), 30 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl create mode 100644 examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl diff --git a/NEWS.md b/NEWS.md index 54fbb90b8fc..5dd911391a6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -38,7 +38,8 @@ for human readability. - Wetting and drying feature and examples for 1D and 2D shallow water equations - Implementation of the polytropic Euler equations in 2D - Implementation of the quasi-1D shallow water equations -- Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` +- Subcell (positivity and local min/max) limiting support for conservative variables + in 2D for `TreeMesh` - AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh` - Added `GradientVariables` type parameter to `AbstractEquationsParabolic` 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 new file mode 100644 index 00000000000..209aa2ae352 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl @@ -0,0 +1,93 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_blast_wave + +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_minmax_variables_cons = ["rho"]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +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 = 6, + n_cells_max = 10_000, + periodicity = false) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.3) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false)) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary 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 new file mode 100644 index 00000000000..c1ba3d96962 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +The Sedov blast wave setup based on Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 +""" +function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + # r0 = 0.5 # = more reasonable setup + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + p0_outer = 1.0e-5 # = true Sedov setup + # p0_outer = 1.0e-3 # = more reasonable setup + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_sedov_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +basis = LobattoLegendreBasis(3) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_minmax_variables_cons = ["rho"]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +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, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.6) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback, + save_solution) +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false)) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 3fea48b30da..282805a0e03 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -39,7 +39,7 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - positivity_variables_cons = [1], + positivity_variables_cons = ["rho"], positivity_correction_factor = 0.5) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, 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 new file mode 100644 index 00000000000..3159a2066ad --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl @@ -0,0 +1,142 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler multicomponent equations + +# 1) Dry Air 2) Helium + 28% Air +equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), + gas_constants = (0.287, 1.578)) + +""" + initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2}) + +A shock-bubble testcase for multicomponent Euler equations +- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman + Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations + [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972) +""" +function initial_condition_shock_bubble(x, t, + equations::CompressibleEulerMulticomponentEquations2D{ + 5, + 2 + }) + # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972 + # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf + # typical domain is rectangular, we change it to a square, as Trixi can only do squares + @unpack gas_constants = equations + + # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving + delta = 0.03 + + # Region I + rho1_1 = delta + rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta + v1_1 = zero(delta) + v2_1 = zero(delta) + p_1 = 101325 + + # Region II + rho1_2 = 1.225 - delta + rho2_2 = delta + v1_2 = zero(delta) + v2_2 = zero(delta) + p_2 = 101325 + + # Region III + rho1_3 = 1.6861 - delta + rho2_3 = delta + v1_3 = -113.5243 + v2_3 = zero(delta) + p_3 = 159060 + + # Set up Region I & II: + inicenter = SVector(zero(delta), zero(delta)) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + if (x[1] > 0.50) + # Set up Region III + rho1 = rho1_3 + rho2 = rho2_3 + v1 = v1_3 + v2 = v2_3 + p = p_3 + elseif (r < 0.25) + # Set up Region I + rho1 = rho1_1 + rho2 = rho2_1 + v1 = v1_1 + v2 = v2_1 + p = p_1 + else + # Set up Region II + rho1 = rho1_2 + rho2 = rho2_2 + v1 = v1_2 + v2 = v2_2 + p = p_2 + end + + return prim2cons(SVector(v1, v2, p, rho1, rho2), equations) +end +initial_condition = initial_condition_shock_bubble + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) + +limiter_idp = SubcellLimiterIDP(equations, basis; + local_minmax_variables_cons = ["rho" * string(i) + for i in eachcomponent(equations)]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.25, -2.225) +coordinates_max = (2.20, 2.225) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 1_000_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.01) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 300 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (Trixi.density,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 600, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false)) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl index 6241ca30938..7856c9bafbd 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl @@ -88,9 +88,8 @@ volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - positivity_variables_cons = [ - (i + 3 for i in eachcomponent(equations))..., - ]) + positivity_variables_cons = ["rho" * string(i) + for i in eachcomponent(equations)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index b2cdff2ab53..fe9ad92467f 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -51,7 +51,7 @@ volume_flux = (flux_derigs_etal, flux_nonconservative_powell_local_symmetric) basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - positivity_variables_cons = [1], + positivity_variables_cons = ["rho"], positivity_correction_factor = 0.5) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index c86f266147c..d7e30ab1621 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -77,15 +77,24 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - (; positivity) = limiter + (; local_minmax, positivity) = limiter (; output_directory) = callback variables = varnames(cons2cons, semi.equations) mkpath(output_directory) open("$output_directory/deviations.txt", "a") do f print(f, "# iter, simu_time") + if local_minmax + for v in limiter.local_minmax_variables_cons + variable_string = string(variables[v]) + print(f, ", " * variable_string * "_min, " * variable_string * "_max") + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end print(f, ", " * string(variables[v]) * "_min") end end @@ -108,15 +117,26 @@ end @inline function finalize_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) - (; positivity) = limiter + (; local_minmax, positivity) = limiter (; idp_bounds_delta) = limiter.cache variables = varnames(cons2cons, semi.equations) println("─"^100) println("Maximum deviation from bounds:") println("─"^100) + if local_minmax + for v in limiter.local_minmax_variables_cons + v_string = string(v) + println("$(variables[v]):") + println("-lower bound: ", idp_bounds_delta[Symbol(v_string, "_min")][2]) + println("-upper bound: ", idp_bounds_delta[Symbol(v_string, "_max")][2]) + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end println(string(variables[v]) * ":\n- positivity: ", idp_bounds_delta[Symbol(string(v), "_min")][2]) end diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 8159becb503..d52eb6edb9e 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,12 +8,35 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - (; positivity) = solver.volume_integral.limiter + (; local_minmax, positivity) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta) = limiter.cache + if local_minmax + for v in limiter.local_minmax_variables_cons + v_string = string(v) + key_min = Symbol(v_string, "_min") + key_max = Symbol(v_string, "_max") + deviation_min = idp_bounds_delta[key_min] + deviation_max = idp_bounds_delta[key_max] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + var = u[v, i, j, element] + deviation_min[1] = max(deviation_min[1], + variable_bounds[key_min][i, j, element] - var) + deviation_max[1] = max(deviation_max[1], + var - variable_bounds[key_max][i, j, element]) + end + deviation_min[2] = max(deviation_min[2], deviation_min[1]) + deviation_max[2] = max(deviation_max[2], deviation_max[1]) + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end key = Symbol(string(v), "_min") deviation = idp_bounds_delta[key] for element in eachelement(solver, cache), j in eachnode(solver), @@ -30,10 +53,19 @@ # 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[Symbol(v_string, "_min")][1], ", ", + idp_bounds_delta[Symbol(v_string, "_max")][1]) + end + end if positivity for v in limiter.positivity_variables_cons - key = Symbol(string(v), "_min") - print(f, ", ", idp_bounds_delta[key][1]) + if v in limiter.local_minmax_variables_cons + continue + end + print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1]) end end println(f) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 78b1b829b06..ba2ad1cd1cd 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -42,6 +42,18 @@ Common choices of the `conversion_function` are [`cons2cons`](@ref) and """ function varnames end +# Return the index of `varname` in `varnames(solution_variables, equations)` if available. +# Otherwise, throw an error. +function get_variable_index(varname, equations; + solution_variables = cons2cons) + index = findfirst(==(varname), varnames(solution_variables, equations)) + if isnothing(index) + throw(ArgumentError("$varname is no valid variable.")) + end + + return index +end + # Add methods to show some information on systems of equations. function Base.show(io::IO, equations::AbstractEquations) # Since this is not performance-critical, we can use `@nospecialize` to reduce latency. @@ -211,8 +223,8 @@ end """ NonConservativeLocal() -Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". -When the argument `nonconservative_type` is of type `NonConservativeLocal`, +Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". +When the argument `nonconservative_type` is of type `NonConservativeLocal`, the function returns the local part of the non-conservative term. """ struct NonConservativeLocal end @@ -220,8 +232,8 @@ struct NonConservativeLocal end """ NonConservativeSymmetric() -Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". -When the argument `nonconservative_type` is of type `NonConservativeSymmetric`, +Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". +When the argument `nonconservative_type` is of type `NonConservativeSymmetric`, the function returns the symmetric part of the non-conservative term. """ struct NonConservativeSymmetric end diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 91ad59b76b6..9e5ebc7f9b5 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -191,6 +191,12 @@ end A subcell limiting volume integral type for DG methods based on subcell blending approaches with a low-order FV method. Used with limiter [`SubcellLimiterIDP`](@ref). +!!! note + Subcell limiting methods are not fully functional on non-conforming meshes. This is + mainly because the implementation assumes that low- and high-order schemes have the same + surface terms, which is not guaranteed for non-conforming meshes. The low-order scheme + with a high-order mortar is not invariant domain preserving. + !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 97843db7743..2fc62f548d2 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -463,4 +463,21 @@ end return nothing end + +""" + get_boundary_outer_state(boundary_condition::BoundaryConditionDirichlet, + cache, t, equations, dg, indices...) +For subcell limiting, the calculation of local bounds for non-periodic domains require the boundary +outer state. This function returns the boundary value at time `t` and for node with spatial +indices `indices`. +""" +@inline function get_boundary_outer_state(boundary_condition::BoundaryConditionDirichlet, + cache, t, equations, dg, indices...) + (; node_coordinates) = cache.elements + + x = get_node_coords(node_coordinates, equations, dg, indices...) + u_outer = boundary_condition.boundary_value_function(x, t, equations) + + return u_outer +end end # @muladd diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 55d402954bf..055e7ce24a4 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -14,12 +14,17 @@ end """ SubcellLimiterIDP(equations::AbstractEquations, basis; - positivity_variables_cons = [], + local_minmax_variables_cons = String[], + positivity_variables_cons = String[], positivity_correction_factor = 0.1) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: -- positivity limiting for conservative variables (`positivity_variables_cons`) +- Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) +- Positivity limiting for conservative variables (`positivity_variables_cons`) + +Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` +and `positivity_variables_cons = ["rho"]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. @@ -41,62 +46,95 @@ The bounds are calculated using the low-order FV solution. The positivity limite This is an experimental feature and may change in future releases. """ struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter + local_minmax::Bool + local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables positivity::Bool positivity_variables_cons::Vector{Int} # Positivity for conservative variables positivity_correction_factor::RealT cache::Cache end -# this method is used when the indicator is constructed as for shock-capturing volume integrals +# this method is used when the limiter is constructed as for shock-capturing volume integrals function SubcellLimiterIDP(equations::AbstractEquations, basis; - positivity_variables_cons = [], + local_minmax_variables_cons = String[], + positivity_variables_cons = String[], positivity_correction_factor = 0.1) + local_minmax = (length(local_minmax_variables_cons) > 0) positivity = (length(positivity_variables_cons) > 0) - bound_keys = Tuple(Symbol(string(v), "_min") for v in positivity_variables_cons) + local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, + equations) + positivity_variables_cons_ = get_variable_index.(positivity_variables_cons, + equations) + + bound_keys = () + if local_minmax + for v in local_minmax_variables_cons_ + v_string = string(v) + bound_keys = (bound_keys..., Symbol(v_string, "_min"), + Symbol(v_string, "_max")) + end + end + for v in positivity_variables_cons_ + if !(v in local_minmax_variables_cons_) + bound_keys = (bound_keys..., Symbol(string(v), "_min")) + end + end cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys) - SubcellLimiterIDP{typeof(positivity_correction_factor), typeof(cache)}(positivity, - positivity_variables_cons, - positivity_correction_factor, - cache) + SubcellLimiterIDP{typeof(positivity_correction_factor), + typeof(cache)}(local_minmax, local_minmax_variables_cons_, + positivity, positivity_variables_cons_, + positivity_correction_factor, cache) end function Base.show(io::IO, limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - @unpack positivity = limiter + (; local_minmax, positivity) = limiter print(io, "SubcellLimiterIDP(") - if !(positivity) + if !(local_minmax || positivity) print(io, "No limiter selected => pure DG method") else print(io, "limiter=(") + local_minmax && print(io, "min/max limiting, ") positivity && print(io, "positivity") print(io, "), ") end + print(io, "Local bounds with FV solution") print(io, ")") end function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - @unpack positivity = limiter + (; local_minmax, positivity) = limiter if get(io, :compact, false) show(io, limiter) else - if !(positivity) + if !(local_minmax || positivity) setup = ["limiter" => "No limiter selected => pure DG method"] else setup = ["limiter" => ""] + if local_minmax + setup = [ + setup..., + "" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)", + ] + end if positivity - string = "positivity with conservative variables $(limiter.positivity_variables_cons)" + string = "positivity for conservative variables $(limiter.positivity_variables_cons)" setup = [setup..., "" => string] setup = [ setup..., "" => " positivity correction factor = $(limiter.positivity_correction_factor)", ] end + setup = [ + setup..., + "Local bounds" => "FV solution", + ] end summary_box(io, "SubcellLimiterIDP", setup) end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 0a72b79ea3f..bc69e55f264 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -30,6 +30,10 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE @unpack alpha = limiter.cache.subcell_limiter_coefficients alpha .= zero(eltype(alpha)) + if limiter.local_minmax + @trixi_timeit timer() "local min/max limiting" idp_local_minmax!(alpha, limiter, + u, t, dt, semi) + end if limiter.positivity @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, semi) end @@ -52,6 +56,173 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE return nothing end +@inline function calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + # Calc bounds inside elements + @threaded for element in eachelement(dg, cache) + var_min[:, :, element] .= typemax(eltype(var_min)) + var_max[:, :, element] .= typemin(eltype(var_max)) + # Calculate bounds at Gauss-Lobatto nodes using u + for j in eachnode(dg), i in eachnode(dg) + var = u[variable, i, j, element] + var_min[i, j, element] = min(var_min[i, j, element], var) + var_max[i, j, element] = max(var_max[i, j, element], var) + + if i > 1 + var_min[i - 1, j, element] = min(var_min[i - 1, j, element], var) + var_max[i - 1, j, element] = max(var_max[i - 1, j, element], var) + end + if i < nnodes(dg) + var_min[i + 1, j, element] = min(var_min[i + 1, j, element], var) + var_max[i + 1, j, element] = max(var_max[i + 1, j, element], var) + end + if j > 1 + var_min[i, j - 1, element] = min(var_min[i, j - 1, element], var) + var_max[i, j - 1, element] = max(var_max[i, j - 1, element], var) + end + if j < nnodes(dg) + var_min[i, j + 1, element] = min(var_min[i, j + 1, element], var) + var_max[i, j + 1, element] = max(var_max[i, j + 1, element], var) + end + end + end + + # Values at element boundary + calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, mesh) +end + +@inline function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, + mesh::TreeMesh2D) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + # Calc bounds at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get neighboring element ids + left = cache.interfaces.neighbor_ids[1, interface] + right = cache.interfaces.neighbor_ids[2, interface] + + orientation = cache.interfaces.orientations[interface] + + for i in eachnode(dg) + index_left = (nnodes(dg), i) + index_right = (1, i) + if orientation == 2 + index_left = reverse(index_left) + index_right = reverse(index_right) + end + var_left = u[variable, index_left..., left] + var_right = u[variable, index_right..., right] + + var_min[index_right..., right] = min(var_min[index_right..., right], + var_left) + var_max[index_right..., right] = max(var_max[index_right..., right], + var_left) + + var_min[index_left..., left] = min(var_min[index_left..., left], var_right) + var_max[index_left..., left] = max(var_max[index_left..., left], var_right) + end + end + + # Calc bounds at physical boundaries + for boundary in eachboundary(dg, cache) + element = cache.boundaries.neighbor_ids[boundary] + orientation = cache.boundaries.orientations[boundary] + neighbor_side = cache.boundaries.neighbor_sides[boundary] + + for i in eachnode(dg) + if neighbor_side == 2 # Element is on the right, boundary on the left + index = (1, i) + boundary_index = 1 + else # Element is on the left, boundary on the right + index = (nnodes(dg), i) + boundary_index = 2 + end + if orientation == 2 + index = reverse(index) + boundary_index += 2 + end + u_outer = get_boundary_outer_state(boundary_conditions[boundary_index], + cache, t, equations, dg, + index..., element) + var_outer = u_outer[variable] + + var_min[index..., element] = min(var_min[index..., element], var_outer) + var_max[index..., element] = max(var_max[index..., element], var_outer) + end + end + + return nothing +end + +@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi) + for variable in limiter.local_minmax_variables_cons + idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) + end + + return nothing +end + +@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + variable_string = string(variable) + var_min = variable_bounds[Symbol(variable_string, "_min")] + var_max = variable_bounds[Symbol(variable_string, "_max")] + calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) + + @threaded for element in eachelement(dg, semi.cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + var = u[variable, i, j, element] + # Real Zalesak type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + + Qp = max(0, (var_max[i, j, element] - var) / dt) + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pp and Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, element] + + Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + + max(0, val_flux2_local) + max(0, val_flux2_local_jp1) + 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 + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qp = abs(Qp) / + (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, element])) + Qm = abs(Qm) / + (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, element])) + + # Calculate alpha at nodes + alpha[i, j, element] = max(alpha[i, j, element], 1 - min(1, Qp, Qm)) + end + end + + return nothing +end + @inline function idp_positivity!(alpha, limiter, u, dt, semi) # Conservative variables for variable in limiter.positivity_variables_cons @@ -79,6 +250,13 @@ end end # Compute bound + if limiter.local_minmax && + variable in limiter.local_minmax_variables_cons && + var_min[i, j, element] >= positivity_correction_factor * var + # Local limiting is more restrictive that positivity limiting + # => Skip positivity limiting for this node + continue + end var_min[i, j, element] = positivity_correction_factor * var # Real one-sided Zalesak-type limiter diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 4a1dc42772d..04a295537a3 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -313,6 +313,34 @@ end end end +@trixi_testset "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), + l2=[ + 0.3517507570120483, + 0.19252291020146015, + 0.19249751956580294, + 0.618717827188004, + ], + linf=[ + 1.6699566795772216, + 1.3608007992899402, + 1.361864507190922, + 2.44022884092527, + ], + tspan=(0.0, 0.5), + initial_refinement_level=4, + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), l2=[ @@ -339,6 +367,34 @@ end end end +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + l2=[ + 0.4328635350273501, + 0.15011135840723572, + 0.15011135840723572, + 0.616129927549474, + ], + linf=[ + 1.6145297181778906, + 0.8614006163026988, + 0.8614006163026972, + 6.450225090647602, + ], + tspan=(0.0, 1.0), + initial_refinement_level=4, + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_sedov_blast_wave.jl (HLLE)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), l2=[ diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 2e808af6473..30d52b37b96 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -75,6 +75,35 @@ end end end +@trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"), + l2=[ + 73.10832638093902, + 1.4599215762968585, + 57176.014861335476, + 0.17812843581838675, + 0.010123079422717837, + ], + linf=[ + 214.50568817511956, + 25.40392579616452, + 152862.41011222568, + 0.564195553101797, + 0.0956331651771212, + ], + initial_refinement_level=3, + tspan=(0.0, 0.001)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_eulermulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index 1d768c0df69..0c5f2bf0a4c 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -414,7 +414,7 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], 0.1, "cache") + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], 0.1, "cache") @test_nowarn show(stdout, limiter_idp) # TODO: TrixiShallowWater: move unit test From 3f3ac9a9903759519992108357194a61117be3df Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 15 Nov 2023 22:55:34 +0100 Subject: [PATCH 3/8] HLLC Flux for Ideal MHD 1D (#1702) * start HLLC MHD 1D Cartesian * Continue HLLC * tests and format * remove ode diff eq and plots * Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Hendrik Ranocha * comments * further 1d simplifications * label variables in hll flux comp * update test vals for updated ode * fmt --------- Co-authored-by: Hendrik Ranocha --- src/equations/ideal_glm_mhd_1d.jl | 142 ++++++++++++++++++++++++++++++ test/test_tree_1d_mhd.jl | 24 +++++ test/test_unit.jl | 1 + 3 files changed, 167 insertions(+) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 85030e8a5ad..a465571989b 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -259,6 +259,148 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat return SVector(f1, f2, f3, f4, f5, f6, f7, f8) end +""" + flux_hllc(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D) + +- Li (2005) +An HLLC Riemann solver for magneto-hydrodynamics +[DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020). +""" +function flux_hllc(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations1D) + # Unpack left and right states + rho_ll, v1_ll, v2_ll, v3_ll, p_ll, B1_ll, B2_ll, B3_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr = cons2prim(u_rr, equations) + + # Conserved variables + rho_v1_ll = u_ll[2] + rho_v2_ll = u_ll[3] + rho_v3_ll = u_ll[4] + + rho_v1_rr = u_rr[2] + rho_v2_rr = u_rr[3] + rho_v3_rr = u_rr[4] + + # Obtain left and right fluxes + f_ll = flux(u_ll, orientation, equations) + f_rr = flux(u_rr, orientation, equations) + + SsL, SsR = min_max_speed_naive(u_ll, u_rr, orientation, equations) + sMu_L = SsL - v1_ll + sMu_R = SsR - v1_rr + if SsL >= 0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + f5 = f_ll[5] + f6 = f_ll[6] + f7 = f_ll[7] + f8 = f_ll[8] + elseif SsR <= 0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + f5 = f_rr[5] + f6 = f_rr[6] + f7 = f_rr[7] + f8 = f_rr[8] + else + # Compute the "HLLC-speed", eq. (14) from paper mentioned above + #= + SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr - B1_ll^2 + B1_rr^2 ) / + (rho_rr * sMu_R - rho_ll * sMu_L) + =# + # Simplification for 1D: B1 is constant + SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr) / + (rho_rr * sMu_R - rho_ll * sMu_L) + + Sdiff = SsR - SsL + + # Compute HLL values for vStar, BStar + # These correspond to eq. (28) and (30) from the referenced paper + # and the classic HLL intermediate state given by (2) + rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff + + v1Star = (SsR * rho_v1_rr - SsL * rho_v1_ll - (f_rr[2] - f_ll[2])) / + (Sdiff * rho_HLL) + v2Star = (SsR * rho_v2_rr - SsL * rho_v2_ll - (f_rr[3] - f_ll[3])) / + (Sdiff * rho_HLL) + v3Star = (SsR * rho_v3_rr - SsL * rho_v3_ll - (f_rr[4] - f_ll[4])) / + (Sdiff * rho_HLL) + + #B1Star = (SsR * B1_rr - SsL * B1_ll - (f_rr[6] - f_ll[6])) / Sdiff + # 1D B1 = constant => B1_ll = B1_rr = B1Star + B1Star = B1_ll + + B2Star = (SsR * B2_rr - SsL * B2_ll - (f_rr[7] - f_ll[7])) / Sdiff + B3Star = (SsR * B3_rr - SsL * B3_ll - (f_rr[8] - f_ll[8])) / Sdiff + if SsL <= SStar + SdiffStar = SsL - SStar + + densStar = rho_ll * sMu_L / SdiffStar # (19) + + mom_1_Star = densStar * SStar # (20) + mom_2_Star = densStar * v2_ll - + (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar # (21) + mom_3_Star = densStar * v3_ll - + (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar # (22) + + #pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2 # (17) + # 1D B1 = constant => B1_ll = B1_rr = B1Star + pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll # (17) + + enerStar = u_ll[5] * sMu_L / SdiffStar + + (pstar * SStar - p_ll * v1_ll - (B1Star * + (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) - + B1_ll * (B1_ll * v1_ll + B2_ll * v2_ll + B3_ll * v3_ll))) / + SdiffStar # (23) + + # Classic HLLC update (32) + f1 = f_ll[1] + SsL * (densStar - u_ll[1]) + f2 = f_ll[2] + SsL * (mom_1_Star - u_ll[2]) + f3 = f_ll[3] + SsL * (mom_2_Star - u_ll[3]) + f4 = f_ll[4] + SsL * (mom_3_Star - u_ll[4]) + f5 = f_ll[5] + SsL * (enerStar - u_ll[5]) + f6 = f_ll[6] + SsL * (B1Star - u_ll[6]) + f7 = f_ll[7] + SsL * (B2Star - u_ll[7]) + f8 = f_ll[8] + SsL * (B3Star - u_ll[8]) + else # SStar <= Ssr + SdiffStar = SsR - SStar + + densStar = rho_rr * sMu_R / SdiffStar # (19) + + mom_1_Star = densStar * SStar # (20) + mom_2_Star = densStar * v2_rr - + (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar # (21) + mom_3_Star = densStar * v3_rr - + (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar # (22) + + #pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2 # (17) + # 1D B1 = constant => B1_ll = B1_rr = B1Star + pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr # (17) + + enerStar = u_rr[5] * sMu_R / SdiffStar + + (pstar * SStar - p_rr * v1_rr - (B1Star * + (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) - + B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) / + SdiffStar # (23) + + # Classic HLLC update (32) + f1 = f_rr[1] + SsR * (densStar - u_rr[1]) + f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2]) + f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3]) + f4 = f_rr[4] + SsR * (mom_3_Star - u_rr[4]) + f5 = f_rr[5] + SsR * (enerStar - u_rr[5]) + f6 = f_rr[6] + SsR * (B1Star - u_rr[6]) + f7 = f_rr[7] + SsR * (B2Star - u_rr[7]) + f8 = f_rr[8] + SsR * (B3Star - u_rr[8]) + end + end + return SVector(f1, f2, f3, f4, f5, f6, f7, f8) +end + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D) diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index b34bdf0660c..447572eee88 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -206,6 +206,30 @@ end end end +@trixi_testset "elixir_mhd_torrilhon_shock_tube.jl (HLLC)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"), + surface_flux=flux_hllc, + l2=[ + 0.4574266553239646, 0.4794143154876439, 0.3407079689595056, + 0.44797768430829343, 0.9206916204424165, + 1.3216517820475193e-16, 0.2889748702415378, + 0.25529778018020927, + ], + linf=[ + 1.217943947570543, 0.8868438459815245, 0.878215340656725, + 0.9710882819266371, 1.6742759645320984, + 2.220446049250313e-16, 0.704710220504591, 0.6562122176458641, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index 0c5f2bf0a4c..aa7840ec7aa 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -958,6 +958,7 @@ end for u in u_values @test flux_hlle(u, u, 1, equations) ≈ flux(u, 1, equations) + @test flux_hllc(u, u, 1, equations) ≈ flux(u, 1, equations) end equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =# From 8bdd0cded6f4252e090e4349a604a0fd23edb0d6 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 16 Nov 2023 15:45:44 +0100 Subject: [PATCH 4/8] CompatHelper: bump compat for Aqua to 0.8 for package test, (keep existing compat) (#1740) * CompatHelper: bump compat for Aqua to 0.8 for package test, (keep existing compat) * only v0.8 of Aqua.jl * update Aqua.jl keyword * format --------- Co-authored-by: CompatHelper Julia Co-authored-by: Hendrik Ranocha Co-authored-by: Hendrik Ranocha --- test/Project.toml | 2 +- test/test_aqua.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index fcb96b9e48f..ecae0ac0900 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,7 +12,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -Aqua = "0.7" +Aqua = "0.8" CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10" Downloads = "1" ForwardDiff = "0.10" diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 9f57791406f..93457caba28 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -11,8 +11,8 @@ include("test_trixi.jl") ambiguities = false, # exceptions necessary for adding a new method `StartUpDG.estimate_h` # in src/solvers/dgmulti/sbp.jl - piracy = (treat_as_own = [Trixi.StartUpDG.RefElemData, - Trixi.StartUpDG.MeshData],)) + piracies = (treat_as_own = [Trixi.StartUpDG.RefElemData, + Trixi.StartUpDG.MeshData],)) end end #module From c8e1f99dc52666cee2d5d8000d881e7f3dac60e0 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 17 Nov 2023 08:36:37 +0100 Subject: [PATCH 5/8] HLL Davis wave speeds & Cartesian Flux Winters for polytropic Euler (#1733) * hll_davis for polytropic euler * flux winters for tree and tests * fmt * typo * remove copied comment * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * orientation_or_normal_direction * typo --------- Co-authored-by: Andrew Winters --- .../elixir_eulerpolytropic_convergence.jl | 63 +++++++++++ src/equations/polytropic_euler_2d.jl | 102 +++++++++++++++++- test/test_structured_2d.jl | 23 ++++ test/test_tree_2d_eulerpolytropic.jl | 35 ++++++ test/test_tree_2d_part2.jl | 3 + test/test_unit.jl | 48 +++++++++ 6 files changed, 272 insertions(+), 2 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl create mode 100644 test/test_tree_2d_eulerpolytropic.jl diff --git a/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl new file mode 100644 index 00000000000..95ef38eb701 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 1.4 +kappa = 0.5 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +initial_condition = initial_condition_convergence_test + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + periodicity = true, + n_cells_max = 30_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.1) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index d4902bbafb2..f5d2f7b0bad 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat return prim2cons(SVector(rho, v1, v2), equations) end -# Calculate 1D flux for a single point in the normal direction +# Calculate 2D flux for a single point in the normal direction # Note, this directional vector is not normalized @inline function flux(u, normal_direction::AbstractVector, equations::PolytropicEulerEquations2D) @@ -135,8 +135,28 @@ end return SVector(f1, f2, f3) end +# Calculate 2D flux for a single point +@inline function flux(u, orientation::Integer, equations::PolytropicEulerEquations2D) + _, v1, v2 = cons2prim(u, equations) + p = pressure(u, equations) + + rho_v1 = u[2] + rho_v2 = u[3] + + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + end + return SVector(f1, f2, f3) +end + """ - flux_winters_etal(u_ll, u_rr, normal_direction, + flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction, equations::PolytropicEulerEquations2D) Entropy conserving two-point flux for isothermal or polytropic gases. @@ -178,6 +198,37 @@ For details see Section 3.2 of the following reference return SVector(f1, f2, f3) end +@inline function flux_winters_etal(u_ll, u_rr, orientation::Integer, + equations::PolytropicEulerEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + p_ll = equations.kappa * rho_ll^equations.gamma + p_rr = equations.kappa * rho_rr^equations.gamma + + # Compute the necessary mean values + if equations.gamma == 1.0 # isothermal gas + rho_mean = ln_mean(rho_ll, rho_rr) + else # equations.gamma > 1 # polytropic gas + rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) + end + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + + if orientation == 1 # x-direction + f1 = rho_mean * 0.5 * (v1_ll + v1_rr) + f2 = f1 * v1_avg + p_avg + f3 = f1 * v2_avg + else # y-direction + f1 = rho_mean * 0.5 * (v2_ll + v2_rr) + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + end + + return SVector(f1, f2, f3) +end + @inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, equations::PolytropicEulerEquations2D) rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) @@ -196,6 +247,53 @@ end return lambda_min, lambda_max end +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::PolytropicEulerEquations2D) + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + # Pressure for polytropic Euler + p_ll = equations.kappa * rho_ll^equations.gamma + p_rr = equations.kappa * rho_rr^equations.gamma + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + if orientation == 1 # x-direction + λ_min = min(v1_ll - c_ll, v1_rr - c_rr) + λ_max = max(v1_ll + c_ll, v1_rr + c_rr) + else # y-direction + λ_min = min(v2_ll - c_ll, v2_rr - c_rr) + λ_max = max(v2_ll + c_ll, v2_rr + c_rr) + end + + return λ_min, λ_max +end + +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector, + equations::PolytropicEulerEquations2D) + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + # Pressure for polytropic Euler + p_ll = equations.kappa * rho_ll^equations.gamma + p_rr = equations.kappa * rho_rr^equations.gamma + + norm_ = norm(normal_direction) + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # The v_normals are already scaled by the norm + λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr) + λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr) + + return λ_min, λ_max +end + @inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D) rho, v1, v2 = cons2prim(u, equations) c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1)) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index dd2248e10b2..17ba001d5ec 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -598,6 +598,29 @@ end end end +@trixi_testset "elixir_eulerpolytropic_convergence.jl: HLL(Davis)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), + solver=DGSEM(polydeg = 3, + surface_flux = FluxHLL(min_max_speed_davis), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)), + l2=[ + 0.0016689832177644243, 0.0025920263793104445, + 0.003281074494629298, + ], + linf=[ + 0.01099488320190023, 0.013309526619350365, + 0.02008032661117909, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_eulerpolytropic_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"), l2=[ diff --git a/test/test_tree_2d_eulerpolytropic.jl b/test/test_tree_2d_eulerpolytropic.jl new file mode 100644 index 00000000000..545cf7274ff --- /dev/null +++ b/test/test_tree_2d_eulerpolytropic.jl @@ -0,0 +1,35 @@ +module TestExamples2DEulerMulticomponent + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") + +@testset "Polytropic Euler" begin +#! format: noindent + +@trixi_testset "elixir_eulerpolytropic_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulerpolytropic_convergence.jl"), + l2=[ + 0.0016689832177626373, 0.0025920263793094526, + 0.003281074494626679, + ], + linf=[ + 0.010994883201896677, 0.013309526619350365, + 0.02008032661117376, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end +end + +end # module diff --git a/test/test_tree_2d_part2.jl b/test/test_tree_2d_part2.jl index d8e86d14f18..622f12109ff 100644 --- a/test/test_tree_2d_part2.jl +++ b/test/test_tree_2d_part2.jl @@ -26,6 +26,9 @@ isdir(outdir) && rm(outdir, recursive = true) # Compressible Euler Multicomponent include("test_tree_2d_eulermulti.jl") + # Compressible Polytropic Euler + include("test_tree_2d_eulerpolytropic.jl") + # Compressible Euler coupled with acoustic perturbation equations include("test_tree_2d_euleracoustics.jl") diff --git a/test/test_unit.jl b/test/test_unit.jl index aa7840ec7aa..20368481607 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -782,6 +782,54 @@ end end end +@timed_testset "Consistency check for HLL flux with Davis wave speed estimates: Polytropic CEE" begin + flux_hll = FluxHLL(min_max_speed_davis) + + gamma = 1.4 + kappa = 0.5 # Scaling factor for the pressure. + equations = PolytropicEulerEquations2D(gamma, kappa) + u = SVector(1.1, -0.5, 2.34) + + orientations = [1, 2] + for orientation in orientations + @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_hll(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + +@timed_testset "Consistency check for Winters flux: Polytropic CEE" begin + for gamma in [1.4, 1.0, 5 / 3] + kappa = 0.5 # Scaling factor for the pressure. + equations = PolytropicEulerEquations2D(gamma, kappa) + u = SVector(1.1, -0.5, 2.34) + + orientations = [1, 2] + for orientation in orientations + @test flux_winters_etal(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_winters_etal(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + end +end + @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin flux_hll = FluxHLL(min_max_speed_davis) From 4fecc09b927bc516258513a673882e6f8b446353 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Nov 2023 09:47:31 +0100 Subject: [PATCH 6/8] set version to v0.6.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index da618dc78c1..e324de5fc77 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.1-pre" +version = "0.6.1" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 09bdf97df69f55075f853842a97fb2cc2faf703d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Nov 2023 09:47:49 +0100 Subject: [PATCH 7/8] set development version to v0.6.2-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e324de5fc77..17c02a2adac 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.1" +version = "0.6.2-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 1635d318b2bf68f9da2bd9392e8ef782c7e562c6 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 18 Nov 2023 16:20:32 +0100 Subject: [PATCH 8/8] Throw error in `trixi_include` when assignment is not found (#1737) * Throw error in `trixi_include` when assignment is not found * Add tests for `trixi_include` * Reformat * Update test/test_unit.jl Co-authored-by: Michael Schlottke-Lakemper * Add more comments * Fix convergence tests and `trixi_include` with `skip_coverage` * Fix tests trying to override non-existing assignments * Add `maxiters` to restart elixirs and exclude for coverage tests * Reformat * Reformat * Add comment and more specific error in `convergence_test` * Update src/auxiliary/special_elixirs.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/auxiliary/special_elixirs.jl Co-authored-by: Michael Schlottke-Lakemper * Fix tests * Use `@test_nowarn_mod` * Fix tests * Fix tests (again) --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- .../elixir_advection_restart.jl | 2 +- .../elixir_advection_restart.jl | 2 +- .../elixir_advection_restart.jl | 2 +- .../elixir_advection_restart.jl | 2 +- .../tree_2d_dgsem/elixir_advection_restart.jl | 2 +- .../tree_3d_dgsem/elixir_advection_restart.jl | 2 +- .../elixir_euler_restart.jl | 2 +- src/auxiliary/special_elixirs.jl | 47 +++++++-- test/test_mpi_p4est_2d.jl | 5 +- test/test_mpi_p4est_3d.jl | 5 +- test/test_p4est_2d.jl | 9 +- test/test_p4est_3d.jl | 9 +- test/test_special_elixirs.jl | 13 +-- test/test_structured_2d.jl | 15 ++- test/test_structured_3d.jl | 5 +- test/test_threaded.jl | 5 +- test/test_tree_3d_advection.jl | 5 +- test/test_trixi.jl | 2 +- test/test_unit.jl | 97 +++++++++++++++++++ test/test_unstructured_2d.jl | 5 +- test/test_visualization.jl | 5 +- 21 files changed, 196 insertions(+), 45 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/examples/p4est_2d_dgsem/elixir_advection_restart.jl index 4f43e122ab3..0f573714c1f 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_restart.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_restart.jl @@ -31,7 +31,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/examples/p4est_3d_dgsem/elixir_advection_restart.jl b/examples/p4est_3d_dgsem/elixir_advection_restart.jl index cd97d69d692..b3dead42399 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_restart.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_restart.jl @@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/examples/structured_2d_dgsem/elixir_advection_restart.jl b/examples/structured_2d_dgsem/elixir_advection_restart.jl index 19863faae8d..0accbdba702 100644 --- a/examples/structured_2d_dgsem/elixir_advection_restart.jl +++ b/examples/structured_2d_dgsem/elixir_advection_restart.jl @@ -30,7 +30,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/examples/structured_3d_dgsem/elixir_advection_restart.jl b/examples/structured_3d_dgsem/elixir_advection_restart.jl index e81ad5d6430..e516d794df8 100644 --- a/examples/structured_3d_dgsem/elixir_advection_restart.jl +++ b/examples/structured_3d_dgsem/elixir_advection_restart.jl @@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/examples/tree_2d_dgsem/elixir_advection_restart.jl b/examples/tree_2d_dgsem/elixir_advection_restart.jl index 770629bb15e..e0d1003f524 100644 --- a/examples/tree_2d_dgsem/elixir_advection_restart.jl +++ b/examples/tree_2d_dgsem/elixir_advection_restart.jl @@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, alg, dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks; ode_default_options()...) + callback = callbacks, maxiters = 100_000; ode_default_options()...) # Load saved context for adaptive time integrator if integrator.opts.adaptive diff --git a/examples/tree_3d_dgsem/elixir_advection_restart.jl b/examples/tree_3d_dgsem/elixir_advection_restart.jl index f81e013fc0a..a53f4ebf5b1 100644 --- a/examples/tree_3d_dgsem/elixir_advection_restart.jl +++ b/examples/tree_3d_dgsem/elixir_advection_restart.jl @@ -27,7 +27,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl index 4d6af65b8a4..c1908691902 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl @@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, callback = callbacks, maxiters = 100_000); # Get the last time index and work with that. load_timestep!(integrator, restart_filename) diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl index 25bca8939ce..5fdd9aea0c5 100644 --- a/src/auxiliary/special_elixirs.jl +++ b/src/auxiliary/special_elixirs.jl @@ -20,7 +20,7 @@ providing examples with sensible default values for users. Before replacing assignments in `elixir`, the keyword argument `maxiters` is inserted into calls to `solve` and `Trixi.solve` with it's default value used in the SciML ecosystem -for ODEs, see the "Miscellaneous" section of the +for ODEs, see the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). # Examples @@ -36,6 +36,16 @@ julia> redirect_stdout(devnull) do ``` """ function trixi_include(mod::Module, elixir::AbstractString; kwargs...) + # Check that all kwargs exist as assignments + code = read(elixir, String) + expr = Meta.parse("begin \n$code \nend") + expr = insert_maxiters(expr) + + for (key, val) in kwargs + # This will throw an error when `key` is not found + find_assignment(expr, key) + end + # Print information on potential wait time only in non-parallel case if !mpi_isparallel() @info "You just called `trixi_include`. Julia may now compile the code, please be patient." @@ -243,6 +253,7 @@ end function find_assignment(expr, destination) # declare result to be able to assign to it in the closure local result + found = false # find explicit and keyword assignments walkexpr(expr) do x @@ -250,12 +261,17 @@ function find_assignment(expr, destination) if (x.head === Symbol("=") || x.head === :kw) && x.args[1] === Symbol(destination) result = x.args[2] + found = true # dump(x) end end return x end + if !found + throw(ArgumentError("assignment `$destination` not found in expression")) + end + result end @@ -274,17 +290,28 @@ function extract_initial_resolution(elixir, kwargs) return initial_refinement_level end catch e - if isa(e, UndefVarError) - # get cells_per_dimension from the elixir - cells_per_dimension = eval(find_assignment(expr, :cells_per_dimension)) - - if haskey(kwargs, :cells_per_dimension) - return kwargs[:cells_per_dimension] - else - return cells_per_dimension + # If `initial_refinement_level` is not found, we will get an `ArgumentError` + if isa(e, ArgumentError) + try + # get cells_per_dimension from the elixir + cells_per_dimension = eval(find_assignment(expr, :cells_per_dimension)) + + if haskey(kwargs, :cells_per_dimension) + return kwargs[:cells_per_dimension] + else + return cells_per_dimension + end + catch e2 + # If `cells_per_dimension` is not found either + if isa(e2, ArgumentError) + throw(ArgumentError("`convergence_test` requires the elixir to define " * + "`initial_refinement_level` or `cells_per_dimension`")) + else + rethrow() + end end else - throw(e) + rethrow() end end end diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl index 1edbce8f6c8..da90537fcfd 100644 --- a/test/test_mpi_p4est_2d.jl +++ b/test/test_mpi_p4est_2d.jl @@ -69,7 +69,10 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[4.507575525876275e-6], - linf=[6.21489667023134e-5]) + linf=[6.21489667023134e-5], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) end @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl index 8082930b3b4..75f43650082 100644 --- a/test/test_mpi_p4est_3d.jl +++ b/test/test_mpi_p4est_3d.jl @@ -63,7 +63,10 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[0.002590388934758452], - linf=[0.01840757696885409]) + linf=[0.01840757696885409], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) end @trixi_testset "elixir_advection_cubed_sphere.jl" begin diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 07c6d02bbcd..db34aecc168 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -94,7 +94,10 @@ end @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[4.507575525876275e-6], - linf=[6.21489667023134e-5]) + linf=[6.21489667023134e-5], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -216,8 +219,8 @@ end ], surface_flux=flux_hlle, tspan=(0.0, 0.3)) - # 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] diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 346c61c7448..f2467f30204 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -110,7 +110,10 @@ end @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[0.002590388934758452], - linf=[0.01840757696885409]) + linf=[0.01840757696885409], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -305,8 +308,8 @@ end ], tspan=(0.0, 0.3), surface_flux=flux_hlle) - # 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] diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 85671002ba6..ba670a6025e 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -92,8 +92,7 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", # the convergence test logic @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", - "elixir_advection_basic.jl"), 2, - tspan = (0.0, 0.01)) + "elixir_advection_basic.jl"), 2) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"), 2, @@ -101,12 +100,10 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", tspan = (0.0, 0.1)) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", - "elixir_advection_basic.jl"), 2, - tspan = (0.0, 0.01)) + "elixir_advection_basic.jl"), 2) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", - "elixir_advection_coupled.jl"), 2, - tspan = (0.0, 0.01)) + "elixir_advection_coupled.jl"), 2) @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_extended.jl"), 2, @@ -345,8 +342,8 @@ end end @timed_testset "elixir_euler_ad.jl" begin - @test_trixi_include(joinpath(examples_dir(), "special_elixirs", - "elixir_euler_ad.jl")) + @test_nowarn_mod trixi_include(joinpath(examples_dir(), "special_elixirs", + "elixir_euler_ad.jl")) end end end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 17ba001d5ec..96202e00f58 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -194,7 +194,10 @@ end @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[4.219208035582454e-6], - linf=[3.438434404412494e-5]) + linf=[3.438434404412494e-5], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -211,7 +214,10 @@ end linf=[0.0015194252169410394], rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS) elixir_file="elixir_advection_waving_flag.jl", - restart_file="restart_000021.h5") + restart_file="restart_000021.h5", + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -227,7 +233,10 @@ end l2=[7.841217436552029e-15], linf=[1.0857981180834031e-13], elixir_file="elixir_advection_free_stream.jl", - restart_file="restart_000036.h5") + restart_file="restart_000036.h5", + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_structured_3d.jl b/test/test_structured_3d.jl index 0213e1a9813..a52c459d6be 100644 --- a/test/test_structured_3d.jl +++ b/test/test_structured_3d.jl @@ -62,7 +62,10 @@ end @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[0.0025903889347585777], - linf=[0.018407576968841655]) + linf=[0.018407576968841655], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_threaded.jl b/test/test_threaded.jl index 478c90b476a..dbbcbf4c7ce 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -208,7 +208,10 @@ end linf=[0.0015194252169410394], rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS) elixir_file="elixir_advection_waving_flag.jl", - restart_file="restart_000021.h5") + restart_file="restart_000021.h5", + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_tree_3d_advection.jl b/test/test_tree_3d_advection.jl index 56278629417..ae53a2df52f 100644 --- a/test/test_tree_3d_advection.jl +++ b/test/test_tree_3d_advection.jl @@ -27,7 +27,10 @@ end @trixi_testset "elixir_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2=[0.00016017848135651983], - linf=[0.0014175368788298393]) + linf=[0.0014175368788298393], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_trixi.jl b/test/test_trixi.jl index 245efbc0175..cebe2164ae6 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -33,7 +33,7 @@ macro test_trixi_include(elixir, args...) local kwargs = Pair{Symbol, Any}[] for arg in args if (arg.head == :(=) && - !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override)) + !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override, :skip_coverage)) && !(coverage && arg.args[1] in keys(coverage_override))) push!(kwargs, Pair(arg.args...)) end diff --git a/test/test_unit.jl b/test/test_unit.jl index 20368481607..7f774e11d3b 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1374,6 +1374,103 @@ end @test mesh.boundary_faces[:entire_boundary] == [1, 2] end + +@testset "trixi_include" begin + @trixi_testset "Basic" begin + example = """ + x = 4 + """ + + filename = tempname() + try + open(filename, "w") do file + write(file, example) + end + + # Use `@trixi_testset`, which wraps code in a temporary module, and call + # `trixi_include` with `@__MODULE__` in order to isolate this test. + @test_warn "You just called" trixi_include(@__MODULE__, filename) + @test @isdefined x + @test x == 4 + + @test_warn "You just called" trixi_include(@__MODULE__, filename, x = 7) + @test x == 7 + + @test_throws "assignment `y` not found in expression" trixi_include(@__MODULE__, + filename, + y = 3) + finally + rm(filename, force = true) + end + end + + @trixi_testset "With `solve` Without `maxiters`" begin + # `trixi_include` assumes this to be the `solve` function of OrdinaryDiffEq, + # and therefore tries to insert the kwarg `maxiters`, which will fail here. + example = """ + solve() = 0 + x = solve() + """ + + filename = tempname() + try + open(filename, "w") do file + write(file, example) + end + + # Use `@trixi_testset`, which wraps code in a temporary module, and call + # `trixi_include` with `@__MODULE__` in order to isolate this test. + @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__, + filename) + + @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__, + filename, + maxiters = 3) + finally + rm(filename, force = true) + end + end + + @trixi_testset "With `solve` with `maxiters`" begin + # We need another example file that we include with `Base.include` first, in order to + # define the `solve` method without `trixi_include` trying to insert `maxiters` kwargs. + # Then, we can test that `trixi_include` inserts the kwarg in the `solve()` call. + example1 = """ + solve(; maxiters=0) = maxiters + """ + + example2 = """ + x = solve() + """ + + filename1 = tempname() + filename2 = tempname() + try + open(filename1, "w") do file + write(file, example1) + end + open(filename2, "w") do file + write(file, example2) + end + + # Use `@trixi_testset`, which wraps code in a temporary module, and call + # `Base.include` and `trixi_include` with `@__MODULE__` in order to isolate this test. + Base.include(@__MODULE__, filename1) + @test_warn "You just called" trixi_include(@__MODULE__, filename2) + @test @isdefined x + # This is the default `maxiters` inserted by `trixi_include` + @test x == 10^5 + + @test_warn "You just called" trixi_include(@__MODULE__, filename2, + maxiters = 7) + # Test that `maxiters` got overwritten + @test x == 7 + finally + rm(filename1, force = true) + rm(filename2, force = true) + end + end +end end end #module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index d4416ac5b6a..5341d86a7d1 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -129,7 +129,10 @@ end 0.005243995459478956, 0.004685630332338153, 0.01750217718347713, - ]) + ], + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 100_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_visualization.jl b/test/test_visualization.jl index 48164a70fb3..6444dc91d5d 100644 --- a/test/test_visualization.jl +++ b/test/test_visualization.jl @@ -243,7 +243,6 @@ end @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl"), - tspan = (0, 0.1), analysis_callback = Trixi.TrivialCallback()) @test adapt_to_mesh_level(sol, 5) isa Tuple @@ -259,7 +258,6 @@ end @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_basic.jl"), - tspan = (0, 0.1), analysis_callback = Trixi.TrivialCallback(), initial_refinement_level = 1) @test PlotData2D(sol) isa Trixi.PlotData2DCartesian @@ -288,8 +286,7 @@ end @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "structured_3d_dgsem", - "elixir_advection_basic.jl"), - tspan = (0, 0.1)) + "elixir_advection_basic.jl")) @testset "1D plot from 3D solution and general mesh" begin @testset "Create 1D plot as slice" begin