PERK4 standalone

examples/structured_2d_dgsem/elixir_euler_vortex_perk4.jl

@@ -0,0 +1,115 @@
+
+using OrdinaryDiffEq # Required for `CallbackSet`
+using Trixi
+
+# Ratio of specific heats
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hllc)
+
+"""
+    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
+
+The classical isentropic vortex test case as presented in Section 5.1 of
+
+- Brian Vermeire (2019).
+  Paired Explicit Runge-Kutta Schemes for Stiff Systems of Equations
+  [DOI:10.1016/j.jcp.2019.05.014](https://doi.org/10.1016/j.jcp.2019.05.014)
+  https://spectrum.library.concordia.ca/id/eprint/985444/1/Paired-explicit-Runge-Kutta-schemes-for-stiff-sy_2019_Journal-of-Computation.pdf
+"""
+function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
+    # Evaluate error after full domain traversion
+    if t == T_end
+        t = 0
+    end
+
+    # initial center of the vortex
+    inicenter = SVector(0.0, 0.0)
+    # strength of the vortex
+    S = 13.5
+    # Radius of vortex
+    R = 1.5
+    # Free-stream Mach 
+    M = 0.4
+    # base flow
+    v1 = 1.0
+    v2 = 1.0
+    vel = SVector(v1, v2)
+
+    center = inicenter + vel * t # advection of center
+    center = x - center          # distance to centerpoint
+    center = SVector(center[2], -center[1])
+    r2 = center[1]^2 + center[2]^2
+
+    f = (1 - r2) / (2 * R^2)
+
+    rho = (1 - (S * M / pi)^2 * (gamma - 1) * exp(2 * f) / 8)^(1 / (gamma - 1))
+
+    du = S / (2 * π * R) * exp(f) # vel. perturbation
+    vel = vel + du * center
+    v1, v2 = vel
+
+    p = rho^gamma / (gamma * M^2)
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+initial_condition = initial_condition_isentropic_vortex
+
+EdgeLength = 20
+
+N_passes = 1
+T_end = EdgeLength * N_passes
+tspan = (0.0, T_end)
+
+coordinates_min = (-EdgeLength / 2, -EdgeLength / 2)
+coordinates_max = (EdgeLength / 2, EdgeLength / 2)
+
+cells_per_dimension = (32, 32)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_integrals = (entropy,))
+
+# Note quite large CFL number
+stepsize_callback = StepsizeCallback(cfl = 9.1)
+
+callbacks = CallbackSet(summary_callback,
+                        stepsize_callback,
+                        analysis_callback)
+
+###############################################################################
+# set up time integration algorithm
+
+num_stages = 19
+
+current_directory = @__DIR__
+coefficient_file = "a_" * string(num_stages) * ".txt"
+
+# Download the optimized PERK4 coefficients
+Trixi.download("https://gist.githubusercontent.com/DanielDoehring/84f266ff61f0a69a0127cec64056275e/raw/1a66adbe1b425d33daf502311ecbdd4b191b89cc/a_19.txt",
+               joinpath(current_directory, coefficient_file))
+
+ode_algorithm = Trixi.PairedExplicitRK4(num_stages, current_directory)
+
+###############################################################################
+# run the simulation
+
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 42.0, # Not used
+                  save_everystep = false, callback = callbacks);
+
+# Clean up the downloaded file                  
+rm(joinpath(current_directory, coefficient_file))
+
+summary_callback()

examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl

@@ -0,0 +1,65 @@
+
+using OrdinaryDiffEq # Required for `CallbackSet`
+using Trixi
+
+# Convex and ECOS are imported because they are used for finding the optimal time step and optimal 
+# monomial coefficients in the stability polynomial of P-ERK time integrators.
+using Convex, ECOS
+
+###############################################################################
+# semidiscretization of the hyperbolic diffusion equations
+
+equations = HyperbolicDiffusionEquations1D()
+
+initial_condition = initial_condition_poisson_nonperiodic
+
+boundary_conditions = boundary_condition_poisson_nonperiodic
+
+solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000,
+                periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions,
+                                    source_terms = source_terms_poisson_nonperiodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+resid_tol = 5.0e-12
+steady_state_callback = SteadyStateCallback(abstol = resid_tol, reltol = 0.0)
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+# Construct fourth order paired explicit Runge-Kutta method with 11 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.PairedExplicitRK4(11, tspan, semi)
+
+cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
+stepsize_callback = StepsizeCallback(cfl = 0.9 * cfl_number)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+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);
+summary_callback() # print the timer summary

ext/TrixiConvexECOSExt.jl

@@ -15,7 +15,8 @@ end
 using LinearAlgebra: eigvals
 
 # Use functions that are to be extended and additional symbols that are not exported
-using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd
+using Trixi: Trixi, undo_normalization!, undo_normalization_PERK4!,
+             bisect_stability_polynomial, bisect_stability_polynomial_PERK4, @muladd
 
 # By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
 # Since these FMAs can increase the performance of many numerical algorithms,
@@ -28,10 +29,14 @@ using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, @muladd
 # relative to consistency order.
 function Trixi.undo_normalization!(gamma_opt, consistency_order, num_stage_evals)
     for k in (consistency_order + 1):num_stage_evals
-        gamma_opt[k - consistency_order] = gamma_opt[k - consistency_order] /
-                                           factorial(k)
+        gamma_opt[k - consistency_order] /= factorial(k)
+    end
+end
+
+function Trixi.undo_normalization_PERK4!(gamma_opt, num_stage_evals)
+    for k in 1:(num_stage_evals - 5)
+        gamma_opt[k] /= factorial(k + 4)
     end
-    return gamma_opt
 end
 
 # Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency
@@ -63,6 +68,44 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals,
     return maximum(abs(pnoms))
 end
 
+# Specialized form of the stability polynomials for fourth-order PERK schemes.
+function stability_polynomials_PERK4!(pnoms, num_stage_evals,
+                                      normalized_powered_eigvals,
+                                      gamma, dt, cS3)
+    num_eig_vals = length(pnoms)
+
+    # Constants arising from the particular form of Butcher tableau chosen for the 4th order PERK methods
+    k1 = 0.001055026310046423 / cS3
+    k2 = 0.03726406530405851 / cS3
+    # Note: `cS3` = c_{S-3} is in principle free, while the other abscissae are fixed to 1.0
+
+    # Initialize with zero'th order (z^0) coefficient
+    for i in 1:num_eig_vals
+        pnoms[i] = 1.0
+    end
+
+    # First `consistency_order` = 4 terms of the exponential
+    for k in 1:4
+        for i in 1:num_eig_vals
+            pnoms[i] += dt^k * normalized_powered_eigvals[i, k]
+        end
+    end
+
+    # "Fixed" term due to choice of the PERK4 Butcher tableau
+    # Required to un-do the normalization of the eigenvalues here
+    pnoms += k1 * dt^5 * normalized_powered_eigvals[:, 5] * factorial(5)
+
+    # Contribution from free coefficients
+    for k in 1:(num_stage_evals - 5)
+        pnoms += (k2 * dt^(k + 4) * normalized_powered_eigvals[:, k + 4] * gamma[k] +
+                  k1 * dt^(k + 5) * normalized_powered_eigvals[:, k + 5] * gamma[k] *
+                  (k + 5))
+    end
+
+    # For optimization only the maximum is relevant
+    return maximum(abs(pnoms))
+end
+
 #=
 The following structures and methods provide a simplified implementation to 
 discover optimal stability polynomial for a given set of `eig_vals`
@@ -158,6 +201,93 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
         gamma_opt = [gamma_opt]
     end
 
+    undo_normalization!(gamma_opt, consistency_order, num_stage_evals)
+
+    return gamma_opt, dt
+end
+
+# Specialized routine for PERK4.
+# For details, see Section 4 in 
+# - D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024).
+# Fourth-Order Paired-Explicit Runge-Kutta Methods
+# [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470)
+function Trixi.bisect_stability_polynomial_PERK4(num_eig_vals,
+                                                 num_stage_evals,
+                                                 dtmax, dteps, eig_vals, cS3;
+                                                 verbose = false)
+    dtmin = 0.0
+    dt = -1.0
+    abs_p = -1.0
+
+    # Construct stability polynomial for each eigenvalue
+    pnoms = ones(Complex{Float64}, num_eig_vals, 1)
+
+    # Init datastructure for monomial coefficients
+    gamma = Variable(num_stage_evals - 5)
+
+    normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals)
+
+    for j in 1:num_stage_evals
+        fac_j = factorial(j)
+        for i in 1:num_eig_vals
+            normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j
+        end
+    end
+
+    if verbose
+        println("Start optimization of stability polynomial \n")
+    end
+
+    # Bisection on timestep
+    while dtmax - dtmin > dteps
+        dt = 0.5 * (dtmax + dtmin)
+
+        # Use last optimal values for gamma in (potentially) next iteration
+        problem = minimize(stability_polynomials_PERK4!(pnoms,
+                                                        num_stage_evals,
+                                                        normalized_powered_eigvals,
+                                                        gamma, dt, cS3))
+
+        solve!(problem,
+               # Parameters taken from default values for EiCOS
+               MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99,
+                                           "delta" => 2e-7,
+                                           "feastol" => 1e-9,
+                                           "abstol" => 1e-9,
+                                           "reltol" => 1e-9,
+                                           "feastol_inacc" => 1e-4,
+                                           "abstol_inacc" => 5e-5,
+                                           "reltol_inacc" => 5e-5,
+                                           "nitref" => 9,
+                                           "maxit" => 100,
+                                           "verbose" => 3); silent = true)
+
+        abs_p = problem.optval
+
+        if abs_p < 1
+            dtmin = dt
+        else
+            dtmax = dt
+        end
+    end
+
+    if verbose
+        println("Concluded stability polynomial optimization \n")
+    end
+
+    gamma_opt = []
+
+    if num_stage_evals > 5
+        gamma_opt = evaluate(gamma)
+
+        # Catch case S = 6 (only one opt. variable)
+        if isa(gamma_opt, Number)
+            gamma_opt = [gamma_opt]
+        end
+
+        undo_normalization_PERK4!(gamma_opt, num_stage_evals)
+    end
+
     return gamma_opt, dt
 end
 end # @muladd