diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index fedf097ec7..50e064068c 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.24.3 + uses: crate-ci/typos@v1.25.0 diff --git a/NEWS.md b/NEWS.md index 9371cafa07..37711d8e57 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.9 from v0.8.x + +#### Changed + +- We removed the first argument `semi` corresponding to a `Semidiscretization` from the + `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]). + The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069]) + ## Changes in the v0.8 lifecycle #### Changed @@ -11,6 +19,8 @@ for human readability. - The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product of the Jacobian and the conserved variables instead of the conserved variables only to make AMR fully conservative ([#2028]). This may change AMR results slightly. +- Subcell (IDP) limiting is now officially supported and not marked as experimental + anymore (see `VolumeIntegralSubcellLimiting`). ## Changes when updating to v0.8 from v0.7.x diff --git a/README.md b/README.md index 5cda16118e..30e9da243d 100644 --- a/README.md +++ b/README.md @@ -37,6 +37,7 @@ installation and postprocessing procedures. Its features include: * Kinetic energy-preserving and entropy-stable methods based on flux differencing * Entropy-stable shock capturing * Positivity-preserving limiting + * Subcell invariant domain-preserving (IDP) limiting * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl) * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/) * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods) diff --git a/docs/src/index.md b/docs/src/index.md index 869caaed85..0e4749dde3 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -30,6 +30,7 @@ installation and postprocessing procedures. Its features include: * Kinetic energy-preserving and entropy-stable methods based on flux differencing * Entropy-stable shock capturing * Positivity-preserving limiting + * Subcell invariant domain-preserving (IDP) limiting * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl) * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/) * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 45b750728c..638aff1b10 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -85,11 +85,11 @@ analysis_interval = 2000 l_inf = 1.0 # Length of airfoil force_boundary_names = (:AirfoilBottom, :AirfoilTop) -drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, +drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names, DragCoefficientPressure(aoa(), rho_inf(), u_inf(equations), l_inf)) -lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, +lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names, LiftCoefficientPressure(aoa(), rho_inf(), u_inf(equations), l_inf)) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 21df724da4..4fb62d358d 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -89,11 +89,11 @@ rho_inf = 1.4 u_inf = 0.38 l_inf = 1.0 # Diameter of circle -drag_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,), +drag_coefficient = AnalysisSurfaceIntegral((:x_neg,), DragCoefficientPressure(aoa, rho_inf, u_inf, l_inf)) -lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,), +lift_coefficient = AnalysisSurfaceIntegral((:x_neg,), LiftCoefficientPressure(aoa, rho_inf, u_inf, l_inf)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 40baef6ef9..43ce7ac1ed 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -120,23 +120,23 @@ summary_callback = SummaryCallback() analysis_interval = 2000 force_boundary_names = (:AirfoilBottom, :AirfoilTop) -drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, +drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names, DragCoefficientPressure(aoa(), rho_inf(), u_inf(equations), l_inf())) -lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, +lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names, LiftCoefficientPressure(aoa(), rho_inf(), u_inf(equations), l_inf())) -drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, +drag_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names, DragCoefficientShearStress(aoa(), rho_inf(), u_inf(equations), l_inf())) -lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, +lift_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names, LiftCoefficientShearStress(aoa(), rho_inf(), u_inf(equations), diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl index 7461198fbb..22e1d1f3d8 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl +++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl @@ -91,7 +91,7 @@ semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_con # combining both semidiscretizations for Euler + self-gravity parameters = ParametersEulerGravity(background_density = 1.5e7, # aka rho0 gravitational_constant = 6.674e-8, # aka G - cfl = 1.6, + cfl = 0.8, # value as used in the paper resid_tol = 1.0e-4, n_iterations_max = 1000, timestep_gravity = timestep_gravity_carpenter_kennedy_erk54_2N!) @@ -105,7 +105,7 @@ ode = semidiscretize(semi, tspan); summary_callback = SummaryCallback() -stepsize_callback = StepsizeCallback(cfl = 1.0) +stepsize_callback = StepsizeCallback(cfl = 0.5) # value as used in the paper save_solution = SaveSolutionCallback(interval = 10, save_initial_solution = true, diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl new file mode 100644 index 0000000000..2b9602ab4a --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl @@ -0,0 +1,84 @@ + +using Convex, ECOS +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = 1.0 +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = -1.0 # minimum coordinate +coordinates_max = 1.0 # maximum coordinate + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 20.0 +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(alive_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.1, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 4, + med_level = 5, med_threshold = 0.1, + max_level = 6, max_threshold = 0.6) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +# Construct second order paired explicit Runge-Kutta method with 6 stages for given simulation setup. +# Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used +# in calculating the polynomial coefficients in the ODE algorithm. +ode_algorithm = Trixi.PairedExplicitRK2(6, tspan, semi) + +# For Paired Explicit Runge-Kutta methods, we receive an optimized timestep for a given reference semidiscretization. +# To allow for e.g. adaptivity, we reverse-engineer the corresponding CFL number to make it available during the simulation. +cfl_number = Trixi.calculate_cfl(ode_algorithm, ode) +stepsize_callback = StepsizeCallback(cfl = cfl_number) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, + alive_callback, + save_solution, + analysis_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation +sol = Trixi.solve(ode, ode_algorithm, + dt = 1.0, # Manual time step value, will be overwritten by the stepsize_callback when it is specified. + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl index d32b9d6f1f..843d8a7aa4 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl @@ -43,10 +43,7 @@ end ############################################################################### # Specify non-periodic boundary conditions -function inflow(x, t, equations::InviscidBurgersEquation1D) - return initial_condition_rarefaction(coordinate_min, t, equations) -end -boundary_condition_inflow = BoundaryConditionDirichlet(inflow) +boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction) function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t, surface_flux_function, diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 1f0b0e7e04..09c6cee9f9 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -43,10 +43,7 @@ end ############################################################################### # Specify non-periodic boundary conditions -function inflow(x, t, equations::InviscidBurgersEquation1D) - return initial_condition_shock(coordinate_min, t, equations) -end -boundary_condition_inflow = BoundaryConditionDirichlet(inflow) +boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock) function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t, surface_flux_function, diff --git a/src/callbacks_stage/subcell_limiter_idp_correction.jl b/src/callbacks_stage/subcell_limiter_idp_correction.jl index 69125ebecd..d234a6a550 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction.jl @@ -23,9 +23,6 @@ called with [`VolumeIntegralSubcellLimiting`](@ref). - Pazner (2020) Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876) - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. """ struct SubcellLimiterIDPCorrection end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index e22d8b14e9..2b5790e9fb 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -20,7 +20,6 @@ For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary name `:Airfoil` in 2D. -- `semi::Semidiscretization`: Passed in to retrieve boundary condition information - `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries where the quantity of interest is computed - `variable::Variable`: Quantity of interest, like lift or drag @@ -29,8 +28,7 @@ struct AnalysisSurfaceIntegral{Variable, NBoundaries} variable::Variable # Quantity of interest, like lift or drag boundary_symbols::NTuple{NBoundaries, Symbol} # Name(s) of the boundary/boundaries - function AnalysisSurfaceIntegral(semi, - boundary_symbols::NTuple{NBoundaries, Symbol}, + function AnalysisSurfaceIntegral(boundary_symbols::NTuple{NBoundaries, Symbol}, variable) where {NBoundaries} return new{typeof(variable), NBoundaries}(variable, boundary_symbols) end diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 07765eeb9a..c3bac82b4d 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -271,9 +271,6 @@ with a low-order FV method. Used with limiter [`SubcellLimiterIDP`](@ref). 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. """ struct VolumeIntegralSubcellLimiting{VolumeFluxDG, VolumeFluxFV, Limiter} <: AbstractVolumeIntegral diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index f560018f46..21330eec45 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -60,9 +60,6 @@ where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) o - Pazner (2020) Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876) - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. """ struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, LimitingOnesidedVariablesNonlinear, Cache} <: diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 285827850c..ffbc325c82 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -19,9 +19,6 @@ The third-order SSP Runge-Kutta method of Shu and Osher. - Shu, Osher (1988) "Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes" (Eq. 2.18) [DOI: 10.1016/0021-9991(88)90177-5](https://doi.org/10.1016/0021-9991(88)90177-5) - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. """ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP numerator_a::SVector{3, Float64} @@ -133,9 +130,6 @@ end The following structures and methods provide the infrastructure for SSP Runge-Kutta methods of type `SimpleAlgorithmSSP`. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. """ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 23a3ceba76..ad385e6df2 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -68,7 +68,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, a_matrix[:, 1] -= A a_matrix[:, 2] = A - return a_matrix, c + return a_matrix, c, dt_opt end # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2 @@ -76,7 +76,6 @@ end function compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs::AbstractString, bS, cS) - # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) for k in 2:num_stages @@ -107,7 +106,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, end @doc raw""" - PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt, bS = 1.0, cS = 0.5) PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization; verbose = false, bS = 1.0, cS = 0.5) @@ -118,6 +117,7 @@ end - `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing monomial coefficients of the stability polynomial of PERK method. The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks. + - `dt_opt` (`Float64`): Optimal time step size for the simulation setup. - `tspan`: Time span of the simulation. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the @@ -144,16 +144,19 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle b1::Float64 bS::Float64 cS::Float64 + dt_opt::Float64 end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + dt_opt, bS = 1.0, cS = 0.5) + # If the user has the monomial coefficients, they also must have the optimal time step a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs, bS, cS) - return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) end # Constructor that calculates the coefficients with polynomial optimizer from a @@ -171,12 +174,12 @@ end function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, bS = 1.0, cS = 0.5) - a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, - eig_vals, tspan, - bS, cS; - verbose) + a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages, + eig_vals, tspan, + bS, cS; + verbose) - return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS) + return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 @@ -232,6 +235,26 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F, k_higher::uType end +""" + calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode) + +This function computes the CFL number once using the initial condition of the problem and the optimal timestep (`dt_opt`) from the ODE algorithm. +""" +function calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode) + t0 = first(ode.tspan) + u_ode = ode.u0 + semi = ode.p + dt_opt = ode_algorithm.dt_opt + + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + cfl_number = dt_opt / max_dt(u, t0, mesh, + have_constant_speed(equations), equations, + solver, cache) + return cfl_number +end + """ add_tstop!(integrator::PairedExplicitRK2Integrator, t) Add a time stop during the time integration process. diff --git a/test/test_paper_self_gravitating_gas_dynamics.jl b/test/test_paper_self_gravitating_gas_dynamics.jl index 63a7a2b6de..1e097efd51 100644 --- a/test/test_paper_self_gravitating_gas_dynamics.jl +++ b/test/test_paper_self_gravitating_gas_dynamics.jl @@ -215,16 +215,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_jeans_instability.jl"), l2=[ - 10733.63378538114, - 13356.780607423452, + 10733.28239179182, + 13356.0533511341, 1.6722844879795038e-6, - 26834.076821148774 + 26833.19833691448 ], linf=[ - 15194.296424901113, - 18881.481685044182, + 15193.794080890715, + 18880.45819785685, 6.809726988008751e-6, - 37972.99700513482 + 37971.74113135785 ], tspan=(0.0, 0.1), atol=4.0e-6) @@ -242,16 +242,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_jeans_instability.jl"), l2=[ - 10734.598193238024, - 13358.217234481384, + 10734.653995035567, + 13357.709699808021, 1.911011743371934e-6, - 26836.487841241516 + 26836.62734552835 ], linf=[ - 15195.661004798487, - 18883.512035906537, + 15195.73481107317, + 18882.799120551972, 7.867948710816926e-6, - 37976.408478975296 + 37976.592992473394 ], tspan=(0.0, 0.1), atol=4.0e-6, # the background field is reatively large, so this corresponds to our usual atol diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 115c5f3c69..f061e2e1c3 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -123,6 +123,25 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 end end + +# Testing the second-order paired explicit Runge-Kutta (PERK) method with the optimal CFL number +@trixi_testset "elixir_advection_perk2_optimal_cfl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2_optimal_cfl.jl"), + l2=[0.0009700887119146429], + linf=[0.00137209242077041]) + # 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) + # 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)) < 8000 + end +end end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index a86327d5ef..42f42693d8 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1671,7 +1671,7 @@ end Trixi.download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt", joinpath(path_coeff_file, "gamma_6.txt")) - ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file) + ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file, 42) # dummy optimal time step (dt_opt plays no role in determining `a_matrix`) @test isapprox(ode_algorithm.a_matrix, [0.12405417889682908 0.07594582110317093