From 472c64b3ccf37acee99fd1f746d0adb7cb4c1178 Mon Sep 17 00:00:00 2001 From: Warisa Date: Sat, 7 Sep 2024 12:03:04 +0200 Subject: [PATCH 01/23] Update stepsize_callback CFL to 3.0 and add calculate_cfl --- .../tree_1d_dgsem/elixir_advection_perk2.jl | 2 +- src/callbacks_step/stepsize.jl | 27 +++++++++++++++++++ .../methods_PERK2.jl | 10 ++++--- 3 files changed, 35 insertions(+), 4 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl index 7db461a079e..add6389b9fc 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl @@ -40,7 +40,7 @@ analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 2.5) +stepsize_callback = StepsizeCallback(cfl = 3.0) alive_callback = AliveCallback(alive_interval = analysis_interval) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 8b5cb958318..e4d87ef1521 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -66,6 +66,10 @@ end semi = integrator.p @unpack cfl_number = stepsize_callback + if isa(integrator.alg, AbstractPairedExplicitRKSingle) + cfl_number = calculate_cfl(u_ode, t, integrator.alg.dt_opt, cfl_number, semi) + end + # Dispatch based on semidiscretization dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number, semi) @@ -90,6 +94,29 @@ function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization) solver, cache) end +# For Paired Explicit Runge-Kutta methods, attempt to use optimal time step obtained from optimizing +# the stability polynomials of the Butcher tableau. If this is not possible, fall back to the CFL condition +# that is either specified by the user in stepsize_callback or the default one. +function calculate_cfl(u_ode, t, dt_opt, default_cfl_number, semi::AbstractSemidiscretization) + 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, t, mesh, + have_constant_speed(equations), equations, + solver, cache) + + println("cfl calculated: ", cfl_number) + + # Ensure that the default CFL number is not exceeded + if cfl_number <= default_cfl_number && cfl_number > 0.0 + println("cfl (using optimal time step): ", cfl_number) + return cfl_number + else + println("cfl (default): ", default_cfl_number) + return default_cfl_number + end +end + # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have # an integrator at this stage but only the ODE, this method will be used there. It's called in 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 23a3ceba76c..ddbf53cdc0b 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 @@ -77,6 +77,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs::AbstractString, bS, cS) + #TODO: If the user has a specific set of monomial coefficients, they must also have already obtained dt_opt + #TODO: where did I get this monomial in unit test.... # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) for k in 2:num_stages @@ -106,6 +108,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, return a_matrix, c end +#TODO: add dt_opt to docstring @doc raw""" PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, bS = 1.0, cS = 0.5) @@ -144,6 +147,7 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle b1::Float64 bS::Float64 cS::Float64 + dt_opt::Float64 end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file @@ -171,12 +175,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, + 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 From 12900690e0d428301158e9bf544a9018283474ac Mon Sep 17 00:00:00 2001 From: Warisa Date: Sun, 8 Sep 2024 22:04:07 +0700 Subject: [PATCH 02/23] delete unnecessary line --- src/callbacks_step/stepsize.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index e4d87ef1521..5324209fb0d 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -105,8 +105,6 @@ function calculate_cfl(u_ode, t, dt_opt, default_cfl_number, semi::AbstractSemid have_constant_speed(equations), equations, solver, cache) - println("cfl calculated: ", cfl_number) - # Ensure that the default CFL number is not exceeded if cfl_number <= default_cfl_number && cfl_number > 0.0 println("cfl (using optimal time step): ", cfl_number) From 4395e54bcf2fd39bdacbeac07ba12e47d37cf8ca Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 9 Sep 2024 16:11:57 +0700 Subject: [PATCH 03/23] update calculation of cfl number --- .../tree_1d_dgsem/elixir_advection_perk2.jl | 3 ++- src/callbacks_step/stepsize.jl | 23 ++++++++----------- 2 files changed, 11 insertions(+), 15 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl index add6389b9fc..af24e6681b0 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl @@ -40,7 +40,8 @@ analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 3.0) +# For PERK schemes, the CFL number is calculated from the optimal time step of the scheme. +stepsize_callback = StepsizeCallback() alive_callback = AliveCallback(alive_interval = analysis_interval) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 5324209fb0d..78429f7c59f 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -64,10 +64,13 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - @unpack cfl_number = stepsize_callback + + # In the case where max_dt varies, calculate all cfl numbers and pick the largest one? if isa(integrator.alg, AbstractPairedExplicitRKSingle) - cfl_number = calculate_cfl(u_ode, t, integrator.alg.dt_opt, cfl_number, semi) + cfl_number = calculate_cfl(u_ode, t, integrator.alg.dt_opt, semi) + else + @unpack cfl_number = stepsize_callback end # Dispatch based on semidiscretization @@ -94,10 +97,9 @@ function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization) solver, cache) end -# For Paired Explicit Runge-Kutta methods, attempt to use optimal time step obtained from optimizing -# the stability polynomials of the Butcher tableau. If this is not possible, fall back to the CFL condition -# that is either specified by the user in stepsize_callback or the default one. -function calculate_cfl(u_ode, t, dt_opt, default_cfl_number, semi::AbstractSemidiscretization) +# For Paired Explicit Runge-Kutta methods, use the CFL number calculated from the optimal timestep of the +# scheme. +function calculate_cfl(u_ode, t, dt_opt, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) @@ -105,14 +107,7 @@ function calculate_cfl(u_ode, t, dt_opt, default_cfl_number, semi::AbstractSemid have_constant_speed(equations), equations, solver, cache) - # Ensure that the default CFL number is not exceeded - if cfl_number <= default_cfl_number && cfl_number > 0.0 - println("cfl (using optimal time step): ", cfl_number) - return cfl_number - else - println("cfl (default): ", default_cfl_number) - return default_cfl_number - end + return cfl_number end # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own From e103667d47107d844e73bf3e788a2cf9ff73a2db Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 9 Sep 2024 16:32:56 +0700 Subject: [PATCH 04/23] update tests --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 8 ++++---- test/test_tree_1d_advection.jl | 5 +++-- test/test_unit.jl | 4 ++-- 3 files changed, 9 insertions(+), 8 deletions(-) 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 ddbf53cdc0b..2073ffaff32 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -60,7 +60,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, eig_vals; verbose) monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) - + num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) @@ -108,9 +108,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, return a_matrix, c end -#TODO: add dt_opt to docstring @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) @@ -121,6 +120,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. - `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 @@ -151,7 +151,7 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file -function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt, bS = 1.0, cS = 0.5) a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs, diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 20586c4f3ba..18fd5985c92 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -82,10 +82,11 @@ end end end +#TODO: replace l2 and linf errors with the values in CI @trixi_testset "elixir_advection_perk2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), - l2=[0.011288030389423475], - linf=[0.01596735472556976]) + l2=[2.11731099e-03], + linf=[3.00261319e-03]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_unit.jl b/test/test_unit.jl index aa66ce2556f..549ee6caa0a 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1670,8 +1670,8 @@ end path_coeff_file = mktempdir() 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 From 844eca1ca186790136143aa891cc740470199b89 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 10 Sep 2024 20:18:48 +0700 Subject: [PATCH 05/23] update cfl number calculation for PERK (put this in the constructor) --- .../tree_1d_dgsem/elixir_advection_perk2.jl | 19 ++++++------ ext/TrixiConvexECOSExt.jl | 2 +- src/callbacks_step/stepsize.jl | 31 ++++++++++++------- 3 files changed, 30 insertions(+), 22 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl index af24e6681b0..d09095977b2 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl @@ -39,10 +39,6 @@ summary_callback = SummaryCallback() analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -# For PERK schemes, the CFL number is calculated from the optimal time step of the scheme. -stepsize_callback = StepsizeCallback() - alive_callback = AliveCallback(alive_interval = analysis_interval) save_solution = SaveSolutionCallback(dt = 0.1, @@ -50,6 +46,15 @@ save_solution = SaveSolutionCallback(dt = 0.1, save_final_solution = true, solution_variables = cons2prim) +# 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) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +# For PERK schemes, the CFL number is calculated from the optimal time step of the scheme. +stepsize_callback = StepsizeCallback(ode, ode_algorithm) + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, alive_callback, @@ -59,12 +64,6 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation - -# 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) - 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); diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 8251fe3eed9..0a5b4c44c5d 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -136,7 +136,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, "reltol_inacc" => 5e-5, "nitref" => 9, "maxit" => 100, - "verbose" => 3); silent = true) + "verbose" => 3); silent_solver = true) abs_p = problem.optval diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 78429f7c59f..5660f835b36 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain", stepsize_callback = cb.affect! setup = [ - "CFL number" => stepsize_callback.cfl_number, + "CFL number" => stepsize_callback.cfl_number ] summary_box(io, "StepsizeCallback", setup) end @@ -47,6 +47,22 @@ function StepsizeCallback(; cfl::Real = 1.0) initialize = initialize!) end +# For Paired-Explicit Runge-Kutta methods, the CFL number is calculated based on the optimal timestep +function StepsizeCallback(ode, ode_algorithm::AbstractPairedExplicitRKSingle) + # TODO: Loop over all the time step and choose the minimum CFL number? from the loop??? + t = first(ode.tspan) + u_ode = ode.u0 + semi = ode.p + dt_opt = ode_algorithm.dt_opt + cfl = calculate_cfl(u_ode, t, dt_opt, semi) + + stepsize_callback = StepsizeCallback(cfl) + + DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! + save_positions = (false, false), + initialize = initialize!) +end + function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, integrator) where {Condition, Affect! <: StepsizeCallback} cb.affect!(integrator) @@ -64,14 +80,7 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - - - # In the case where max_dt varies, calculate all cfl numbers and pick the largest one? - if isa(integrator.alg, AbstractPairedExplicitRKSingle) - cfl_number = calculate_cfl(u_ode, t, integrator.alg.dt_opt, semi) - else - @unpack cfl_number = stepsize_callback - end + @unpack cfl_number = stepsize_callback # Dispatch based on semidiscretization dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number, @@ -104,8 +113,8 @@ function calculate_cfl(u_ode, t, dt_opt, semi::AbstractSemidiscretization) u = wrap_array(u_ode, mesh, equations, solver, cache) cfl_number = dt_opt / max_dt(u, t, mesh, - have_constant_speed(equations), equations, - solver, cache) + have_constant_speed(equations), equations, + solver, cache) return cfl_number end From ceb49bdd08edb8792ea7f6a8402afca7c9feaafd Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 10 Sep 2024 20:20:01 +0700 Subject: [PATCH 06/23] revert an unnecessary change on TrixiConvexECOS --- ext/TrixiConvexECOSExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 0a5b4c44c5d..8251fe3eed9 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -136,7 +136,7 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, "reltol_inacc" => 5e-5, "nitref" => 9, "maxit" => 100, - "verbose" => 3); silent_solver = true) + "verbose" => 3); silent = true) abs_p = problem.optval From 466643cfe6199ee066554ba43b311d17b43c1fbb Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 10 Sep 2024 20:24:40 +0700 Subject: [PATCH 07/23] revert another unnecessary change i made in stepsize.jl --- src/callbacks_step/stepsize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 5660f835b36..a6c74cb5f1e 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain", stepsize_callback = cb.affect! setup = [ - "CFL number" => stepsize_callback.cfl_number + "CFL number" => stepsize_callback.cfl_number, ] summary_box(io, "StepsizeCallback", setup) end From 7c0e26d9edd8e6ba8db1843f2b1e2cabd031c9c0 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 10 Sep 2024 20:27:18 +0700 Subject: [PATCH 08/23] update test values but need to be changed again according to CI workflow --- test/test_tree_1d_advection.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 18fd5985c92..9e311684b5b 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -85,8 +85,8 @@ end #TODO: replace l2 and linf errors with the values in CI @trixi_testset "elixir_advection_perk2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), - l2=[2.11731099e-03], - linf=[3.00261319e-03]) + l2=[1.45247275e-02], + linf=[2.05428436e-02]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -97,6 +97,7 @@ end end end +# TODO: edit these values as well # Testing the second-order paired explicit Runge-Kutta (PERK) method without stepsize callback @trixi_testset "elixir_advection_perk2.jl(fixed time step)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), From 8a70681034b6b8a5bb705bdd7afc44dd6b9a8bae Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 20:40:43 +0700 Subject: [PATCH 09/23] revert changes made in test_tree_1d_advaction.jl since we shouldn't alter the existing example --- test/test_tree_1d_advection.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 9e311684b5b..20586c4f3ba 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -82,11 +82,10 @@ end end end -#TODO: replace l2 and linf errors with the values in CI @trixi_testset "elixir_advection_perk2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), - l2=[1.45247275e-02], - linf=[2.05428436e-02]) + l2=[0.011288030389423475], + linf=[0.01596735472556976]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -97,7 +96,6 @@ end end end -# TODO: edit these values as well # Testing the second-order paired explicit Runge-Kutta (PERK) method without stepsize callback @trixi_testset "elixir_advection_perk2.jl(fixed time step)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2.jl"), From b4c1d478cc86cd768a96595a6ae07e37226154f7 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 20:43:00 +0700 Subject: [PATCH 10/23] revert changes made in stepsize.jl --- src/callbacks_step/stepsize.jl | 31 +------------------------------ 1 file changed, 1 insertion(+), 30 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index a6c74cb5f1e..0a95cb35ad3 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain", stepsize_callback = cb.affect! setup = [ - "CFL number" => stepsize_callback.cfl_number, + "CFL number" => stepsize_callback.cfl_number ] summary_box(io, "StepsizeCallback", setup) end @@ -47,22 +47,6 @@ function StepsizeCallback(; cfl::Real = 1.0) initialize = initialize!) end -# For Paired-Explicit Runge-Kutta methods, the CFL number is calculated based on the optimal timestep -function StepsizeCallback(ode, ode_algorithm::AbstractPairedExplicitRKSingle) - # TODO: Loop over all the time step and choose the minimum CFL number? from the loop??? - t = first(ode.tspan) - u_ode = ode.u0 - semi = ode.p - dt_opt = ode_algorithm.dt_opt - cfl = calculate_cfl(u_ode, t, dt_opt, semi) - - stepsize_callback = StepsizeCallback(cfl) - - DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! - save_positions = (false, false), - initialize = initialize!) -end - function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, integrator) where {Condition, Affect! <: StepsizeCallback} cb.affect!(integrator) @@ -106,19 +90,6 @@ function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization) solver, cache) end -# For Paired Explicit Runge-Kutta methods, use the CFL number calculated from the optimal timestep of the -# scheme. -function calculate_cfl(u_ode, t, dt_opt, semi::AbstractSemidiscretization) - 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, t, mesh, - have_constant_speed(equations), equations, - solver, cache) - - return cfl_number -end - # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have # an integrator at this stage but only the ODE, this method will be used there. It's called in From d994e343a74d468d67dc68fae6c75aac65dacfac Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 20:46:44 +0700 Subject: [PATCH 11/23] bring back the old example --- .../tree_1d_dgsem/elixir_advection_perk2.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl index d09095977b2..7db461a079e 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl @@ -39,6 +39,9 @@ summary_callback = SummaryCallback() analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 2.5) + alive_callback = AliveCallback(alive_interval = analysis_interval) save_solution = SaveSolutionCallback(dt = 0.1, @@ -46,15 +49,6 @@ save_solution = SaveSolutionCallback(dt = 0.1, save_final_solution = true, solution_variables = cons2prim) -# 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) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -# For PERK schemes, the CFL number is calculated from the optimal time step of the scheme. -stepsize_callback = StepsizeCallback(ode, ode_algorithm) - # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, alive_callback, @@ -64,6 +58,12 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation + +# 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) + 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); From fffe9db1a7b317a83966e559125180856d1c0150 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 21:17:18 +0700 Subject: [PATCH 12/23] add new example and create a function that simply return cfl number calculated from dt_opt +fmt --- .../elixir_advection_perk2_optimal_cfl.jl | 73 +++++++++++++++++++ .../methods_PERK2.jl | 31 ++++++-- test/test_unit.jl | 12 +-- 3 files changed, 105 insertions(+), 11 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl 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 00000000000..691c7e038e6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl @@ -0,0 +1,73 @@ + +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) + +# 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, + 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/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 2073ffaff32..9ddb27006e5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -60,7 +60,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, eig_vals; verbose) monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages) - + num_monomial_coeffs = length(monomial_coeffs) @assert num_monomial_coeffs == coeffs_max A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) @@ -151,7 +151,8 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle end # struct PairedExplicitRK2 # Constructor that reads the coefficients from a file -function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt, +function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, + dt_opt, bS = 1.0, cS = 0.5) a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs, @@ -176,9 +177,9 @@ function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64}; verbose = false, bS = 1.0, cS = 0.5) a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages, - eig_vals, tspan, - bS, cS; - verbose) + eig_vals, tspan, + bS, cS; + verbose) return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt) end @@ -236,6 +237,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_unit.jl b/test/test_unit.jl index 549ee6caa0a..34d84d94b52 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1288,7 +1288,7 @@ end 0.5011914484393387, 0.8829127712445113, 0.43024132987932817, - 0.7560616633050348, + 0.7560616633050348 ] equations = CompressibleEulerEquations2D(1.4) @@ -1446,7 +1446,7 @@ end SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, - flux_hllc, flux_chandrashekar, + flux_hllc, flux_chandrashekar ] for f_std in fluxes @@ -1471,7 +1471,7 @@ end SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, - flux_hllc, flux_chandrashekar, + flux_hllc, flux_chandrashekar ] for f_std in fluxes @@ -1510,7 +1510,7 @@ end flux_central, flux_hindenlang_gassner, FluxHLL(min_max_speed_davis), - flux_hlle, + flux_hlle ] for f_std in fluxes @@ -1537,7 +1537,7 @@ end flux_central, flux_hindenlang_gassner, FluxHLL(min_max_speed_davis), - flux_hlle, + flux_hlle ] for f_std in fluxes @@ -1670,7 +1670,7 @@ end path_coeff_file = mktempdir() 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, 42) # dummy optimal time step (dt_opt plays no role in determining a_matrix) @test isapprox(ode_algorithm.a_matrix, From 593cbfbb144235a760d3df2d880b883852690259 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 21:19:41 +0700 Subject: [PATCH 13/23] revert unnecessary change from formatting I made in stepsize.jl --- src/callbacks_step/stepsize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 0a95cb35ad3..8b5cb958318 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain", stepsize_callback = cb.affect! setup = [ - "CFL number" => stepsize_callback.cfl_number + "CFL number" => stepsize_callback.cfl_number, ] summary_box(io, "StepsizeCallback", setup) end From 8478bf806c074700b8a8d00adb82e21261f767bb Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 21:22:54 +0700 Subject: [PATCH 14/23] revert unnecessary fmt changes --- test/test_unit.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 34d84d94b52..8575e6f32cf 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1288,7 +1288,7 @@ end 0.5011914484393387, 0.8829127712445113, 0.43024132987932817, - 0.7560616633050348 + 0.7560616633050348, ] equations = CompressibleEulerEquations2D(1.4) @@ -1446,7 +1446,7 @@ end SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, - flux_hllc, flux_chandrashekar + flux_hllc, flux_chandrashekar, ] for f_std in fluxes @@ -1510,7 +1510,7 @@ end flux_central, flux_hindenlang_gassner, FluxHLL(min_max_speed_davis), - flux_hlle + flux_hlle, ] for f_std in fluxes @@ -1537,7 +1537,7 @@ end flux_central, flux_hindenlang_gassner, FluxHLL(min_max_speed_davis), - flux_hlle + flux_hlle, ] for f_std in fluxes From a98f50e93281fcac9671707b4dad2be8706676a1 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 21:24:01 +0700 Subject: [PATCH 15/23] revert another change in test_unit.jl --- test/test_unit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 8575e6f32cf..d77bae5aa19 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1471,7 +1471,7 @@ end SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, - flux_hllc, flux_chandrashekar + flux_hllc, flux_chandrashekar, ] for f_std in fluxes From f78703f192be7dcd17c97bef0aa3f5b0a20f3336 Mon Sep 17 00:00:00 2001 From: Warisa Date: Wed, 11 Sep 2024 22:28:31 +0700 Subject: [PATCH 16/23] correct a constructor + add and delete some comments. --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) 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 9ddb27006e5..e8153160a70 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -76,9 +76,6 @@ end function compute_PairedExplicitRK2_butcher_tableau(num_stages, base_path_monomial_coeffs::AbstractString, bS, cS) - - #TODO: If the user has a specific set of monomial coefficients, they must also have already obtained dt_opt - #TODO: where did I get this monomial in unit test.... # c Vector form Butcher Tableau (defines timestep per stage) c = zeros(num_stages) for k in 2:num_stages @@ -154,11 +151,12 @@ end # struct PairedExplicitRK2 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 From 303173ed0f2b9f96ff62a63716cbe79d1cccad59 Mon Sep 17 00:00:00 2001 From: Warisa Date: Thu, 12 Sep 2024 16:25:46 +0700 Subject: [PATCH 17/23] add a test for optimal cfl number of PERK2 --- test/test_tree_1d_advection.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 20586c4f3ba..d76cce6d7c0 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -115,6 +115,21 @@ 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=[1.45247275e-02], + linf=[2.05428436e-02]) #TODO: fix these values to match the ones in CI workflow + # 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)) < 8000 + end +end end end # module From 49c9873b7dc756c4da5967ea9fa87b53399b9c39 Mon Sep 17 00:00:00 2001 From: Warisa Date: Fri, 13 Sep 2024 21:49:29 +0700 Subject: [PATCH 18/23] correct the values of test to math thosein CI workflow --- test/test_tree_1d_advection.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index d76cce6d7c0..37f5c6bf05f 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -119,8 +119,8 @@ 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=[1.45247275e-02], - linf=[2.05428436e-02]) #TODO: fix these values to match the ones in CI workflow + l2=[0.014524726892608971], + linf=[0.0205428426617561]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From ef12c85de212263908aa4fdce79d0e6ef8e0dce7 Mon Sep 17 00:00:00 2001 From: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> Date: Tue, 17 Sep 2024 15:14:48 +0700 Subject: [PATCH 19/23] Update test/test_unit.jl Co-authored-by: Daniel Doehring --- test/test_unit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 34d84d94b52..5831122ffe2 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, 42) # dummy optimal time step (dt_opt plays no role in determining a_matrix) + 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 From 2bc740349f696c486dfc0996d39fef290fa5ea02 Mon Sep 17 00:00:00 2001 From: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> Date: Tue, 17 Sep 2024 15:15:07 +0700 Subject: [PATCH 20/23] Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring --- .../paired_explicit_runge_kutta/methods_PERK2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 e8153160a70..ad385e6df24 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -117,7 +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. + - `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 From 15570b585d93aea46bbc98fff820c896f76a8019 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 17 Sep 2024 15:34:54 +0700 Subject: [PATCH 21/23] use amr with the current example --- .../elixir_advection_perk2_optimal_cfl.jl | 13 ++++++++++++- test/test_tree_1d_advection.jl | 2 +- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl index 691c7e038e6..c96ed89826f 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl @@ -46,6 +46,16 @@ save_solution = SaveSolutionCallback(dt = 0.1, 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. @@ -60,7 +70,8 @@ stepsize_callback = StepsizeCallback(cfl = cfl_number) callbacks = CallbackSet(summary_callback, alive_callback, save_solution, - analysis_callback, + analysis_callback, + amr_callback, stepsize_callback) ############################################################################### diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 37f5c6bf05f..f3f793274e5 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -120,7 +120,7 @@ end @trixi_testset "elixir_advection_perk2_optimal_cfl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2_optimal_cfl.jl"), l2=[0.014524726892608971], - linf=[0.0205428426617561]) + linf=[0.0205428426617561]) #TODO: update this value to match the new example # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 4ed6c7428b584c86a164bddacaa637b09c33c933 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 17 Sep 2024 18:10:55 +0700 Subject: [PATCH 22/23] fix test values + fmt --- examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl | 2 +- test/test_tree_1d_advection.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl index c96ed89826f..2b9602ab4a4 100644 --- a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl @@ -70,7 +70,7 @@ stepsize_callback = StepsizeCallback(cfl = cfl_number) callbacks = CallbackSet(summary_callback, alive_callback, save_solution, - analysis_callback, + analysis_callback, amr_callback, stepsize_callback) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 81071513ee9..bf44b7e072d 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -127,8 +127,8 @@ 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.014524726892608971], - linf=[0.0205428426617561]) #TODO: update this value to match the new example + l2=[0.0009700887119146429], + linf=[0.00137209242077041]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 9ab373d7e1f35a4cf272c0eb121f0963a1b48471 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 17 Sep 2024 14:50:15 +0200 Subject: [PATCH 23/23] Update test/test_tree_1d_advection.jl --- test/test_tree_1d_advection.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index bf44b7e072d..f061e2e1c30 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -135,6 +135,10 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # 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