Adaptive MIRK BVP solvers - GSoC final report
Qingyu Qu
  2023-08-22 11:40:36 2023-08-22 11:40   2023-08-24 11:41:36      1.3k Words  8 Mins  

In this Google Summer of Code 2023 project, I developed defect control adaptivity according to MIRKDC FORTRAN program for BoundaryValueDiffEq.jl and the related paper: Runge-Kutta Software with Defect Control for Boundary Value ODEs. After implementing the basic structure of defect control adaptivity, I focused on working on improving the user experience of BoundaryValueDiffEq.jl, and helped maintain OrdinaryDiffEq.jl.

It was a really fruitful summer, I would like to thank my mentors, Chris Rackauckas, Avik Pal and Yingbo Ma, it was their guidance that helped me complete my GSoC project.

What was done

To complete GSoC project, I broke down the whole project into three parts:

Defect Control Adaptivity in MIRK Methods

The detailed defect control adaptivity process can be simplified as: With the discrete solution we obtained from nonlinear solving, we first estimate the defect on a specific time interval. To estimate the defect on a subinterval, we need to evaluate the interpolation weights at sample points and . After we get the interpolation weights, we need to set up the interpolations before defect estimation, basically, the interpolation setup flow is to prepare the extra stages for interpolation construction. After we get everything we need, we can proceed to construct the interpolant. We add the discrete solution, RK method stages(by-product during we evaluate in the collocation process) and extra stages to construct the interpolant at two sample points. Then we can start to evaluate the defect during the whole period using the interpolant.

To evaluate the defect on each subinterval, we simply evaluate and . As for the defect2 we have and . Then the bigger estimated defect is the one we wanted.

After we get the defect on the whole time interval, we can continue our computing, if the defect norm is bigger than the threshold(here, the threshold is predefined as 10%), it means the solution is not acceptable, we set the return code as Failure and start a mesh redistribution. If the defect norm is acceptable, we continue to compare the defect norm with our absolute tolerance, if the tolerance is satisfied, we get our final solution. But if the tolerance is not been satisfied, we need to construct a new mesh to equidistribute the defect and start the nonlinear solving and defect estimation process all over again.

To sum up, the core of the defect control adaptivity is we estimate the defect, and according to the defect, we decide whether we need to refine our mesh to achieve better accuracy. Just as described in New Interpolants for Asymptotically Correct Defect Control of BVODEs, the goal of defect control is to compute a numerical solution by adaptively choosing a mesh so that the maximum defect over the entire problem domain is bounded by a user-provided tolerance.

To experience the defect control adaptivity in BoundaryValueDiffEq, we can simply set adaptive = true when we solve our BVP in the solve function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
using BoundaryValueDiffEq
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g / L) * sin(θ)
end
function bc1!(residual, u, p, t)
    residual[1] = u[1][1] + pi / 2
    residual[2] = u[end][1] - pi / 2
end
L = 1.0
const g = 9.81
tspan = (0.0, pi / 2)
bvp=BVProblem(simplependulum!,bc1!,[0.0, 0.0],tspan)
sol = solve(bvp, MIRK4(), adaptive = true)

Improve the internals of BoundaryValueDiffEq.jl

With several new MIRK methods being developed for BoundaryValueDiffEq.jl, the plan for the near future is to optimize the internals of BoundaryValueDiffEq.jl, for details we can refer to More TODOs

Enhance the BVP ecosystem in SciML

Besides the adaptivity work, I implemented the standard test problems for BVP solvers and worked on the SciMLBase.jl, improving the user experience of BVP solvers, including fixing initial guesses, and specific functions for BVP and many more.

Detailed PRs

BoundaryValueDiffEq.jl

Description Status
Add MIRK6 method Merged
Add defect control adaptivity Merged
Fix stage in collocation Merged
Add MIRK5 method Merged
refactor: Fix MIRKTableau and MIRKInterpTableau Merged
Add MIRK3 method Merged
Update README Merged
Update docs for solvers Merged
Add MIRK2 method Merged
Fix initial guess interface Merged

OrdinaryDiffEq.jl

Description Status
Implement DPRKN5 method Merged
Add ERKN7 method Merged
Add more optional threading Merged
Remove unnecessary branches of Vern7, Vern9 and Tsit5 Merged
Add docs for optional threading Merged
Add DPRKN4 method Merged
Add DPRKN6FM method Merged
Throw error when using partitioned ODE methods with ODEProblem Merged
Add MSRK54 method Merged
Add MSRK6 method Merged
Add ESDIRK436L2SA2 method Merged
Add CONTRIBUTING.md Merged
Add ESDIRK437L2SA method Merged
Add ESDIRK547L2SA2 Merged
Add SIR54 method Merged
Add CG4a method Merged
Add missing methods to docs Merged
Finish all stage_limiter! and step_limiter! Merged
Change exp to exponential! in linear methods Merged

SciMLBase.jl

Description Status
Add isinplace checks for DDEFunction Merged
Add BVPFunction Merged
Throw appropriate error when analytic is not given Merged
Improve initial guess for TwoPointBVProblem Merged
Fix promoted tspan for BVProblem and ODEProblem Merged

DiffEqDocs.jl

Description Status
Fix DiffEqDevDocs broken link Merged
Add docs for DPRKN5 and ERKN7 Merged
Add docs for DPRKN4 and DPRKN6FM Merged
Add docs for MSRK 5, 6 and Stepanov5 methods Merged
Add DiscreteFunctionExpr, DiscreteProblemExpr and docs Merged
Add docs for new ESDIRK methods Merged
Add docs for SIR54, CG4a and Alshina methods Merged
Update BVP tutorial Merged
Add BVP test problems docs Merged
Fix BVProblemLibrary docs Merged
Update BVP docs Merged

DiffEqProblemLibrary.jl

Description Status
Add BVProblemLibrary Merged
Add nonlinear BVP test problems Merged
Add the rest BVP test problems Merged

Other repos

Description Status
Add DiscreteFunction Merged
Add transistor amplifier DAE problem Merged
Revive docs Merged
Replace ADAM with Adam etc Merged
Fix TrustRegion docstring Merged

Weekly Reports

These are the weekly reports that I recorded during GSoC period:

Week Link
Week1: https://erikqqy.github.io/2023/06/04/gsoc-1st-week/
Week2: https://erikqqy.github.io/2023/06/11/gsoc-2nd-week/
Week3: https://erikqqy.github.io/2023/06/18/gsoc-3rd-week/
Week4: https://erikqqy.github.io/2023/06/25/gsoc-4th-week/
Week5: https://erikqqy.github.io/2023/07/02/gsoc-5th-week/
Week6: https://erikqqy.github.io/2023/07/09/gsoc-6th-week/
Week7: https://erikqqy.github.io/2023/07/16/gsoc-7th-week/
Week8: https://erikqqy.github.io/2023/07/23/gsoc-8th-week/
Week9: https://erikqqy.github.io/2023/07/30/gsoc-9th-week/
Week10: https://erikqqy.github.io/2023/08/06/gsoc-10th-week/
Week11: https://erikqqy.github.io/2023/08/13/gsoc-11st-week/
Week12: https://erikqqy.github.io/2023/08/20/gsoc-12nd-week/

Future Work

I will continue to maintain and improve BoundaryValueDiffEq.jl and other DiffEq packages in the future, help implement more algorithms and make the internal of BoundaryValueDiffEq.jl much more robust.

NeuralBVP

With the recent development of BoundaryValueDiffEq.jl, especially the type stability and performance enhancement and defect control adaptivity, the BVP solving stack in SciML is growing stronger, We can now utilize BoundaryValueDiffEq.jl to help research NeuralBVP and related topics.

Right now, there are only MIRK method and shooting method in BoundaryValueDiffEq.jl, while FORTRAN has many other solvers like COLDAE to solve BVP-DAE, MUSL and MUSN to perform multiple-shooting methods on BVP, we should add those methods to BoundaryValueDiffEq.jl as well.

Benchmarks on existing solvers

With the recent development in BoundaryValueDiffEq.jl, we need a comprehensive comparison of current solvers in Julia, FORTRAN and MATLAB.

Optimize the current MIRK Methods and Adaptivity

Although for now, the BVP solving is much better than before, there are still a lot of things we can push forward. The first is the new algorithms for solving BVP-related stuff, including but not limited to BVP-DAE solving, multiple shooting etc.

There is still a long way to go to establish a comprehensive and competitive BVP solvers library, but I believe with our joint effort, the BoundaryValueDiffEq.jl and SciML ecosystem will grow more and more powerful!

Acknowledgement

It was a great opportunity to be part of an inspiring community, many thanks to my mentors Chris Rackauckas, Avik Pal and Yingbo Ma. And finally thanks GSoC community and NumFOCUS community for their support.

Adaptive MIRK BVP solvers - GSoC final report
2023/08/22/gsoc-final-report/
Author
Qingyu Qu
Published
2023-08-22 11:40
License
  1. 1. What was done
    1. 1.1. Defect Control Adaptivity in MIRK Methods
    2. 1.2. Improve the internals of BoundaryValueDiffEq.jl
    3. 1.3. Enhance the BVP ecosystem in SciML
  2. 2. Detailed PRs
    1. 2.1. BoundaryValueDiffEq.jl
    2. 2.2. OrdinaryDiffEq.jl
    3. 2.3. SciMLBase.jl
    4. 2.4. DiffEqDocs.jl
    5. 2.5. DiffEqProblemLibrary.jl
    6. 2.6. Other repos
  3. 3. Weekly Reports
  4. 4. Future Work
    1. 4.1. NeuralBVP
    2. 4.2. More BVP-related solvers
    3. 4.3. Benchmarks on existing solvers
    4. 4.4. Optimize the current MIRK Methods and Adaptivity
  5. 5. Acknowledgement