Perk p2 single ext

.github/workflows/Downgrade.yml

@@ -72,7 +72,7 @@ jobs:
       - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-downgrade-compat@v1
         with:
-          skip: LinearAlgebra,Printf,SparseArrays,UUIDs,DiffEqBase
+          skip: LinearAlgebra,Printf,SparseArrays,UUIDs,DiffEqBase,DelimitedFiles
           projects: ., test
       - uses: julia-actions/julia-buildpkg@v1
         env:

NEWS.md

@@ -7,12 +7,13 @@ for human readability.
 ## Changes in the v0.7 lifecycle
 
 #### Added
-- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
-  to 1D and 3D on `TreeMesh` ([#1855], [#1873]).
+- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh` ([#1855], [#1873]).
 - Implementation of 1D Linearized Euler Equations ([#1867]).
 - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients ([#1812]).
 - Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed ([#1835]).
 - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` ([#1792]).
+- New time integrator `PairedExplicitRK2`, implementing the second-order paired explicit Runge-Kutta
+  method with [Convex.jl](https://github.com/jump-dev/Convex.jl) and [ECOS.jl](https://github.com/jump-dev/ECOS.jl) ([#1908])
 
 ## Changes when updating to v0.7 from v0.6.x
 

Project.toml

@@ -7,6 +7,7 @@ version = "0.7.14-pre"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -50,17 +51,23 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
+ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 
 [extensions]
 TrixiMakieExt = "Makie"
+TrixiConvexECOSExt = ["Convex", "ECOS"]
 
 [compat]
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
+Convex = "0.15.4"
 DataStructures = "0.18.15"
+DelimitedFiles = "1"
 DiffEqBase = "6 - 6.143"
 DiffEqCallbacks = "2.25"
 Downloads = "1.6"
+ECOS = "1.1.2"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"
 ForwardDiff = "0.10.24"
@@ -103,3 +110,5 @@ julia = "1.8"
 
 [extras]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
+ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"

examples/tree_1d_dgsem/elixir_advection_perk2.jl

@@ -0,0 +1,66 @@
+
+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)
+
+# 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)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback,
+                        alive_callback,
+                        analysis_callback,
+                        stepsize_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, # 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()

ext/TrixiConvexECOSExt.jl

@@ -0,0 +1,158 @@
+# Package extension for adding Convex-based features to Trixi.jl
+module TrixiConvexECOSExt
+
+# Required for coefficient optimization in P-ERK scheme integrators
+if isdefined(Base, :get_extension)
+    using Convex: MOI, solve!, Variable, minimize, evaluate
+    using ECOS: Optimizer
+else
+    # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
+    using ..Convex: MOI, solve!, Variable, minimize, evaluate
+    using ..ECOS: Optimizer
+end
+
+# Use other necessary libraries
+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
+
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+# Undo normalization of stability polynomial coefficients by index factorial
+# 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)
+    end
+    return gamma_opt
+end
+
+# Compute stability polynomials for paired explicit Runge-Kutta up to specified consistency
+# order, including contributions from free coefficients for higher orders, and
+# return the maximum absolute value
+function stability_polynomials!(pnoms, consistency_order, num_stage_evals,
+                                normalized_powered_eigvals_scaled,
+                                gamma)
+    num_eig_vals = length(pnoms)
+
+    # Initialize with zero'th order (z^0) coefficient
+    for i in 1:num_eig_vals
+        pnoms[i] = 1.0
+    end
+
+    # First `consistency_order` terms of the exponential
+    for k in 1:consistency_order
+        for i in 1:num_eig_vals
+            pnoms[i] += normalized_powered_eigvals_scaled[i, k]
+        end
+    end
+
+    # Contribution from free coefficients
+    for k in (consistency_order + 1):num_stage_evals
+        pnoms += gamma[k - consistency_order] * normalized_powered_eigvals_scaled[:, k]
+    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`
+These are designed for the one-step (i.e., Runge-Kutta methods) integration of initial value ordinary 
+and partial differential equations.
+
+- Ketcheson and Ahmadia (2012).
+Optimal stability polynomials for numerical integration of initial value problems
+[DOI: 10.2140/camcos.2012.7.247](https://doi.org/10.2140/camcos.2012.7.247)
+=#
+
+# Perform bisection to optimize timestep for stability of the polynomial
+function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals,
+                                           num_stage_evals,
+                                           dtmax, dteps, eig_vals;
+                                           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 - consistency_order)
+
+    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
+
+    normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals)
+
+    if verbose
+        println("Start optimization of stability polynomial \n")
+    end
+
+    # Bisection on timestep
+    while dtmax - dtmin > dteps
+        dt = 0.5 * (dtmax + dtmin)
+
+        # Compute stability polynomial for current timestep
+        for k in 1:num_stage_evals
+            dt_k = dt^k
+            for i in 1:num_eig_vals
+                normalized_powered_eigvals_scaled[i, k] = dt_k *
+                                                          normalized_powered_eigvals[i,
+                                                                                     k]
+            end
+        end
+
+        # Use last optimal values for gamma in (potentially) next iteration
+        problem = minimize(stability_polynomials!(pnoms, consistency_order,
+                                                  num_stage_evals,
+                                                  normalized_powered_eigvals_scaled,
+                                                  gamma))
+
+        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_solver = 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
+
+    return evaluate(gamma), dt
+end
+end # @muladd
+
+end # module TrixiConvexECOSExt

src/Trixi.jl

@@ -310,6 +310,14 @@ function __init__()
         end
     end
 
+    @static if !isdefined(Base, :get_extension)
+        @require Convex="f65535da-76fb-5f13-bab9-19810c17039a" begin
+            @require ECOS="e2685f51-7e38-5353-a97d-a921fd2c8199" begin
+                include("../ext/TrixiConvexECOSExt.jl")
+            end
+        end
+    end
+
     # FIXME upstream. This is a hacky workaround for
     #       https://github.com/trixi-framework/Trixi.jl/issues/628
     #       https://github.com/trixi-framework/Trixi.jl/issues/1185