Perk p2 single ext (#1908)

* Templates for PERK p2, p3 * example * constructor changed * Modified so that both of the constructor stay * bring back eps * correct val * add constructor of PERK3 * minor fixes * function and variable name adjustments * spelling fix * Change names and spacing according to style guide * spelling correction * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * snake case * revisit perk single p2 p3 * fmt * fmt * semantic ordering * add literature * Strip code of p = 3 * Add the line to show the error of the elixir * Make adjustments to a test file and delete example of PERK3 * Update Project.toml Co-authored-by: Daniel Doehring <[email protected]> * Update examples/tree_1d_dgsem/elixir_advection_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Add comments in an example of PERK2 Co-authored-by: Daniel Doehring <[email protected]> * Update src/Trixi.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * add tspan as a parameter and adjust the elixir accordingly * add an additional constructor and modify the function compute_PERK2_butcher_tableau * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * update filter function in polynomial optimizer, add literature, and change se_factors to bc_factors * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * change tspan to have (,) instead of[,]. filter_eigenvalues minor adjustments * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/polynomial_optimizer.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Apply suggestions from code review * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * Apply suggestions from code review * use readdlm instead of read_file * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * add DelimFIles * fmt * Unit tests * compat * ecos compat * compat * compat * compat * compat * callbacks * increase allowed allocs * fmt * timer for step callbacks * deps compat * remove del files compat * skip delimitedfiles in downgrade compat * v1 * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * modular imp of the integrator * modularized PolynamialOptimizer * PolynomialOpt modularized * use "fake" extension * make ECOS weakdep * fix name * comment + fmt * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * comment, fix, and make threshold optional * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * edit tspan and alive_interval * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * fmt * fix fmr in test unit * fix fmt in test_unit.jl * some fixes according to code review * add ECOS as a part of enabling TrixiConvexECOSExt.jl compiles * state functions and classes from Convex package used in polynomial optimization * fmt * remove unconditional output in bisection * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * apply suggestion from code review * minor fix * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * add verbose as an optional argument for printing out some outputs during bisection and filtering eigenvalues * remove the file with upper case * add the file with corrected name * minor fix * move constructors outside of struct * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * alter some names from being abbreviated and some fixes * add comments for some functions and move all polynomial optimizaton related functions to TrixiConvexECOSExt.jl * add comments to functions computing butcher tableau * fmt * apply suggestions + fmt * apply suggestion according to code review * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * fix PERK2's name * add short comment regarding PERK's abbreviation * fix export * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Daniel Doehring <[email protected]> * minor corrections * remove exported of PERK functions in extension * set verbose's default value in kwarg * fix path_monomial_coeffs * fmt * add ECOS as dependency in test/Project.toml * bring back polynomial_optimizer so that user can call the constructor with file path can use filter_eigvals without having to load Convex and ECOS * fix values for some tests and add use Convex and ECOS to load function extensions * fix error undefined b1 * attempt to fix error from test by specifying Trixi.entropy * fmt * adjust the value of tests to allign with one in CI and add print command for a test in test_unit.jl * exclude convex warning * update test vals * more warning exclusions * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * minor fix * add what is used from ECOS * spell check * add information about this PR in NEWS.md * Update NEWS.md Co-authored-by: Joshua Lampert <[email protected]> * fix NEWS.md * change Trixi.entropy back to just entropy * Update NEWS.md Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Apply suggestions from code review * Apply suggestions from code review * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl * add some explanations regarding the arguments of the constructors * fmt * exchange "c_end" for "cS" for consistency * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Joshua Lampert <[email protected]> * add more explaination * minor adjustment * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Joshua Lampert <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Joshua Lampert <[email protected]> * Update src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl Co-authored-by: Joshua Lampert <[email protected]> --------- Co-authored-by: Daniel_Doehring <[email protected]> Co-authored-by: Daniel Doehring <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Joshua Lampert <[email protected]>
trixi-framework · May 23, 2024 · c2513e2 · c2513e2
1 parent 887bab9
commit c2513e2
Show file tree

Hide file tree

Showing 14 changed files with 760 additions and 4 deletions.
diff --git a/.github/workflows/Downgrade.yml b/.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:

diff --git a/NEWS.md b/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
 

diff --git a/Project.toml b/Project.toml
@@ -7,6 +7,7 @@ version = "0.7.15-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"
diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/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()
diff --git a/ext/TrixiConvexECOSExt.jl b/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
diff --git a/src/Trixi.jl b/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