Skip to content

Conversation

GodotMisogi
Copy link
Contributor

@GodotMisogi GodotMisogi commented Jul 7, 2025

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

This merges changes to resolve issue #1284 on the benchmarks of handwritten PDEs as ODEs via finite-difference and pseudo-spectral methods.

The package dependencies have been updated to remove unused ones, and the compat entries have been updated for Julia v1.10. The implementations have been tested on local machines (Macbook Pro M1 Pro on macOS Sequoia and a beefier CPU running Ubuntu 22.04 LTS) by running SciMLBenchmarks.weave_file on the relevant .jmd files.

At present, changes have been implemented for the following PDEs in one spatial and one time dimension:

  1. Allen-Cahn equation
  2. Burgers' equation
  3. Kortewreg-de Vries equation
  4. Kuramoto-Sivashinsky equation

The following changes have been implemented to the finite-difference method benchmarks:

  1. The implementation errors in the PDEs have been corrected; some have been updated to reflect their particular characteristics.
  2. Sparse matrices are constructed using SparseArrays.
  3. The MatrixOperator interfaces are used from the updated SciMLOperators for SplitODEProblem.
  4. Sketches of the equations have been added (but require upcoming refinements and updates).
  5. Space-time heatmaps of the solutions instead of line plots.

These are functioning for certain solver algorithm choices (IMEXEuler, CNAB2, CLAF2, ETDRK2, SBDF2)

Commits will be progressively added in this pull (such as pseudo-spectral method updates), including revisions based on comments from the reviewers (@ChrisRackauckas).

@ChrisRackauckas
Copy link
Member

For the psuedospectral, see the comment in #929. The ApproxFun dev mentioned we should change to OrthogonalPolynomials.jl for this part so we should take that recommendation.

@GodotMisogi
Copy link
Contributor Author

GodotMisogi commented Jul 12, 2025

@ChrisRackauckas Could you please run the workflow so I can check the logs for errors?

  1. Allen-Cahn and Burgers' equations via Chebyshev pseudo-spectral methods have been implemented using ClassicalOrthogonalPolynomials.jl. These have been tested with some mass matrix solvers (Rodas3, Rodas4, Rodas5, GRK4T, GRK4A).
  2. The Kuramoto-Sivashinsky and Korteweg-De Vries equations via Fourier pseudo-spectral methods are now working with ApproxFun.jl. These two cases will likely need to be reformulated to utilize Chebyshev interpolants with Neumann boundary conditions, allowing for the use of mass matrix solvers.

@ChrisRackauckas
Copy link
Member

allowing for the use of mass matrix solvers

Why are mass matrices required here? Most of these should resolve to not have any algebraic constraints if the right spectral discretization is chosen. (Singular) Mass matrices (DAEs) would prohibit most of the solvers from being used. Though if that's a requirement for one of them, then so bit and that changes the solvers that should be chosen there (and would be an interesting contrast)

L = 2.0 # Domain length
xs, prob = allen_cahn(N, L)
# @time sol = solve(prob, AutoVern7(RadauIIA5(autodiff=false)); dt=1e-2, abstol=1e-14, reltol=1e-14, adaptive=true)
@time sol = solve(prob, Tsit5(); abstol=1e-14, reltol=1e-14, adaptive=true)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explicit methods usually have issues with convergence due to order loss. I'm not sure changing to an explicit method for these for the reference will work out well. Though you'll see the in the work-precision diagram if you have a problem with the convergence, since all of the other methods will hit a wall and only converge to the error of the reference solution.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch, I commented the wrong line in this commit. Thanks!

reltols = 0.1 .^ (1:3)
multipliers = 0.5 .^ (0:3)
setups = [
Dict(:alg => IMEXEuler(), :dts => 1e-4 * multipliers),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd keep the same sets as above. A set showing all of the exponential integrators, all of the IMEX methods, all of the normal methods, and then a comparison of the best between the families.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Understood, planning to reincorporate the complete tests for the server builds as I can only prototype a few sets with the current local resources. Should these still be fixed time-steps or use adaptive time-stepping?

Dict(:alg => Rodas4(), :dts => 1e-4 * multipliers),
Dict(:alg => Rodas5(), :dts => 1e-4 * multipliers),
Dict(:alg => GRK4T(), :dts => 1e-4 * multipliers),
Dict(:alg => GRK4A(), :dts => 1e-4 * multipliers),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the main ones to add to the mass matrices comparisons would be FBDF and QNDF. Also a DAEProblem formulation testing DASKR.jl and Sundials.IDA. You can benchmark multiple simultaneous problem implementations using the prob_choice setup, see the DAE benchmarks for how this is done.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great to know; thanks for the tips!

@GodotMisogi
Copy link
Contributor Author

GodotMisogi commented Jul 18, 2025

@ChrisRackauckas Please run the workflow so I can check the logs for errors?

  1. The problem descriptions, including equations, for all partial differential equation benchmark tests have been added.
  2. Allen-Cahn and Burgers' equations via finite-difference methods have been implemented using SummationByPartsOperators.jl for conveniently specifying the accuracy of the derivative operator representations.
  3. Allen-Cahn and Burgers' equations via Chebyshev pseudo-spectral methods have been implemented using ClassicalOrthogonalPolynomials.jl. The use of mass matrix solvers has been side-stepped by implementing LU factorization of the (here, non-singular) mass matrix and inverting the linear systems in the vector field implementations provided to SplitODEProblem.
  4. The Kuramoto-Sivashinsky and Korteweg-De Vries equations via finite-difference and Fourier pseudo-spectral methods are implemented using SummationByPartsOperators.jl with periodic boundary conditions.
  5. All tests have included the previous sets of benchmarks comparing implicit-explicit and exponential integrator methods, along with comparisons between the families. Some solvers have been excluded due to excessively long runtimes and some report instabilities in time-integration, which will need further diagnostics.

Comment on lines 229 to 173
setups = [Dict(:alg => KenCarp5()),
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense)),
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
Dict(:alg => ETDRK4(), :dts => 1e-3 * multipliers)]
labels = hcat("KenCarp5 (dense linsolve)", "ARKODE (dense linsolve)", "KenCarp5 (Krylov linsolve)",
"ARKODE (Krylov linsolve)", "ETDRK3 (m=5)", "ETDRK4 (m=5)")
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
print_names=true, names=labels,
numruns=5, error_estimate=:l2,
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="Between family, medium order")
abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [
Dict(:alg => CNAB2(), :dts => 1e-4 * multipliers),
Dict(:alg => CNAB2(linsolve=KrylovJL_GMRES()), :dts => 1e-4 * multipliers),
Dict(:alg => ETDRK2(), :dts => 1e-4 * multipliers),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is missing all of the good high order methods?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I should've specified that the high-tolerance solvers were tested. Currently testing the low-tolerance/high-order methods; will be included in the next few commits.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChrisRackauckas Please run the workflow on the latest commit? I think it's mostly ready.

@GodotMisogi GodotMisogi reopened this Aug 1, 2025
@ChrisRackauckas
Copy link
Member

Wow this is looking pretty good. Very thorough, everything looks correct. I think there's a few little tweaks to maybe do, like adding preconditioners on the Krylov methods, but what I'm going to do is merge it so it can build and it's easier to see it all together. I think it's like 95% or 99% done, and the last little bits I can probably do most of. But let's see what it looks like in the docs. This is a great change, using all of the latest PDE tools and everything!

@ChrisRackauckas ChrisRackauckas merged commit 0a5ccb8 into SciML:master Aug 1, 2025
1 check passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants