example

examples/tree_1d_dgsem/elixir_advection_PERK2.jl

@@ -0,0 +1,66 @@
+
+using OrdinaryDiffEq
+using Trixi
+using Downloads: download
+
+###############################################################################
+# 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 1.0
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# 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_callback = AnalysisCallback(semi, interval = 100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 2.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+MonCoeffsFile = "gamma_6.txt"
+
+download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt",
+MonCoeffsFile)
+
+ode_algorithm = PERK2(6, "./")
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
+
+using Plots
+plot(sol)

src/Trixi.jl

@@ -274,7 +274,7 @@ export trixi_include, examples_dir, get_examples, default_example,
 
 export ode_norm, ode_unstable_check
 
-export PERK, PERK3
+export PERK2, PERK3
 
 export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure
 

src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl

@@ -29,11 +29,6 @@ function ComputePERK2_ButcherTableau(NumStages::Int, BasePathMonCoeffs::Abstract
   for k in 2:NumStages
     c[k] = cEnd * (k - 1)/(NumStages - 1)
   end
-  #=
-  for k in 2:NumStages
-    c[k] = (k - 1)/(2.0*(NumStages - 1))
-  end
-  =#
   println("Timestep-split: "); display(c); println("\n")
   SE_Factors = bS * reverse(c[2:end-1])
 
@@ -86,7 +81,7 @@ mutable struct PERK2 <: PERKSingle
   cEnd::Float64
 
   # Constructor for previously computed A Coeffs
-  function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64, cEnd_::Float64)
+  function PERK2(NumStages_::Int, BasePathMonCoeffs_::AbstractString, bS_::Float64=1.0, cEnd_::Float64=0.5)
 
     newPERK2 = new(NumStages_)
 
@@ -121,7 +116,7 @@ abstract type PERKSingle_Integrator <: PERK_Integrator end
 # This implements the interface components described at
 # https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1
 # which are used in Trixi.
-mutable struct PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK2_IntegratorOptions} <: PERKSingle_Integrator
+mutable struct PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK_IntegratorOptions} <: PERKSingle_Integrator
   u::uType
   du::uType
   u_tmp::uType
@@ -139,7 +134,6 @@ mutable struct PERK2_Integrator{RealT<:Real, uType, Params, Sol, F, Alg, PERK2_I
   k1::uType
   k_higher::uType
   t_stage::RealT
-  du_ode_hyp::uType # TODO: Not best solution since this is not needed for hyperbolic problems
 end
 
 # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771)
@@ -169,7 +163,7 @@ function solve(ode::ODEProblem, alg::PERK2;
   integrator = PERK2_Integrator(u0, du, u_tmp, t0, dt, zero(dt), iter, ode.p,
                               (prob=ode,), ode.f, alg,
                               PERK_IntegratorOptions(callback, ode.tspan; kwargs...), false,
-                              k1, k_higher, t0, du_ode_hyp)
+                              k1, k_higher, t0)
 
   # initialize callbacks
   if callback isa CallbackSet