examples, more constructors

examples/tree_1d_dgsem/elixir_hypdiff_nonperiodic_perk4.jl

@@ -0,0 +1,65 @@
+
+using OrdinaryDiffEq
+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 third order paired explicit Runge-Kutta method with 8 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

src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl

@@ -49,8 +49,6 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
                                                           dtmax,
                                                           dteps,
                                                           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
@@ -112,13 +110,13 @@ end
        In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
     - `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
+    - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the
       equation has been semidiscretized.
     - `verbose` (`Bool`, optional): Verbosity flag, default is false.
-    - `bS` (`Float64`, optional): Value of b in the Butcher tableau at b_s, when 
-      s is the number of stages, default is 1.0.
-    - `cS` (`Float64`, optional): Value of c in the Butcher tableau at c_s, when
-      s is the number of stages, default is 0.5.
+    - `bS` (`Float64`, optional): Value of $b_S$ in the Butcher tableau, where
+      $S$ is the number of stages. Default is 1.0.
+    - `cS` (`Float64`, optional): Value of $c_S$ in the Butcher tableau, where
+      $S$ is the number of stages. Default is 0.5.
 
 The following structures and methods provide a minimal implementation of
 the second-order paired explicit Runge-Kutta (PERK) method
@@ -146,6 +144,7 @@ end # struct PairedExplicitRK2
 function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString,
                            dt_opt = nothing,
                            bS = 1.0, cS = 0.5)
+    @assert num_stages>=2 "PERK2 requires at least two stages"
     # 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,
@@ -159,6 +158,7 @@ end
 function PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization;
                            verbose = false,
                            bS = 1.0, cS = 0.5)
+    @assert num_stages>=2 "PERK2 requires at least two stages"
     eig_vals = eigvals(jacobian_ad_forward(semi))
 
     return PairedExplicitRK2(num_stages, tspan, eig_vals; verbose, bS, cS)
@@ -169,6 +169,7 @@ end
 function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64};
                            verbose = false,
                            bS = 1.0, cS = 0.5)
+    @assert num_stages>=2 "PERK2 requires at least two stages"
     a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages,
                                                                     eig_vals, tspan,
                                                                     bS, cS;

src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl

@@ -48,8 +48,6 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan,
                                                               num_eig_vals, num_stages,
                                                               dtmax, dteps,
                                                               eig_vals; verbose)
-        monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order,
-                                              num_stages)
 
         # Solve the nonlinear system of equations from monomial coefficient and
         # Butcher array abscissae c to find Butcher matrix A
@@ -114,11 +112,11 @@ end
        In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
     - `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
+    - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the
       equation has been semidiscretized.
     - `verbose` (`Bool`, optional): Verbosity flag, default is false.
-    - `cS2` (`Float64`, optional): Value of c in the Butcher tableau at c_{s-2}, when
-      s is the number of stages, default is 1.0f0.
+    - `cS2` (`Float64`, optional): Value of $c_{S-2}$ in the Butcher tableau, where
+      $S$ is the number of stages. Default is 1.0f0.
 
 The following structures and methods provide an implementation of
 the third-order paired explicit Runge-Kutta (P-ERK) method
@@ -147,6 +145,7 @@ end # struct PairedExplicitRK3
 function PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString,
                            dt_opt = nothing;
                            cS2 = 1.0f0)
+    @assert num_stages>=3 "PERK3 requires at least three stages"
     a_matrix, c = compute_PairedExplicitRK3_butcher_tableau(num_stages,
                                                             base_path_a_coeffs;
                                                             cS2)
@@ -157,6 +156,7 @@ end
 # Constructor that computes Butcher matrix A coefficients from a semidiscretization
 function PairedExplicitRK3(num_stages, tspan, semi::AbstractSemidiscretization;
                            verbose = false, cS2 = 1.0f0)
+    @assert num_stages>=3 "PERK3 requires at least three stages"
     eig_vals = eigvals(jacobian_ad_forward(semi))
 
     return PairedExplicitRK3(num_stages, tspan, eig_vals; verbose, cS2)
@@ -165,6 +165,7 @@ end
 # Constructor that calculates the coefficients with polynomial optimizer from a list of eigenvalues
 function PairedExplicitRK3(num_stages, tspan, eig_vals::Vector{ComplexF64};
                            verbose = false, cS2 = 1.0f0)
+    @assert num_stages>=3 "PERK3 requires at least three stages"
     a_matrix, c, dt_opt = compute_PairedExplicitRK3_butcher_tableau(num_stages,
                                                                     tspan,
                                                                     eig_vals;