SDG 16th biweekly report

Review of the past two weeks

In the last two weeks, my work can be summarized as two parts: Refactoring the boundary conditions convention in BoundaryValueDiffEq.jl and make the code more autodiff compatible.

1. This is kind of old problems with BoundaryValueDiffEq.jl, detailed information can be checked in numerous issues: README example should use interpolation instead of indexing, Indexing in boundary conditions, and previous various attempt on resolving this: Boundary conditions should always use solution object. But finally got some time for this issue now: Boundary conditions should always use solution object and Use solution object in all solvers, the basic idea is to change the previous bc convention from indexing based:

function bc!(residual, u, p, t)
    residual[1] = u[:, end÷2][1] + pi / 2
    residual[2] = u[:, end][1] - pi / 2
end

to interpolation based:

function bc!(residual, u, p, t)
    residual[1] = u(pi/4)[1] + pi / 2
    residual[2] = u(pi/2)[1] - pi / 2
end

The previous methods are overly focused on changing the inner value of a solution, which is difficult since there would be dual numbers in the solution, the current workaround is to changing the solution in the loss function, but in the Jacobian function building process, we just built a new one. Interesting that it does accelerate the solving process:

# Current workaround
julia> @benchmark sol = solve(prob, MIRK4(), dt = 0.05)
BenchmarkTools.Trial: 9976 samples with 1 evaluation.
 Range (min … max):  398.833 μs …  21.430 ms  ┊ GC (min … max):  0.00% … 97.59%
 Time  (median):     434.625 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   500.182 μs ± 563.331 μs  ┊ GC (mean ± σ):  12.65% ± 10.75%

  █▄                                                            ▁
  ██▇▅▄▃▁▁▁▃▁▁▁▃▁▃▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▃▆▇ █
  399 μs        Histogram: log(frequency) by time       4.41 ms <

 Memory estimate: 1.08 MiB, allocs estimate: 24342.

# Master branch
julia> @benchmark sol = solve(prob, MIRK4(), dt = 0.05)
BenchmarkTools.Trial: 9520 samples with 1 evaluation.
 Range (min … max):  411.166 μs …  15.173 ms  ┊ GC (min … max):  0.00% … 96.72%
 Time  (median):     459.167 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   524.214 μs ± 598.530 μs  ┊ GC (mean ± σ):  11.43% ±  9.61%

  █▂                                                            ▁
  ██▄▃▃▁▃▄▄▁▃▁▁▁▄▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▁▃▁▁▃▁▁▁▃▁▃▁▁▃▁▁▁▃▆ █
  411 μs        Histogram: log(frequency) by time       5.38 ms <

 Memory estimate: 1.03 MiB, allocs estimate: 23762.

2. While implementing changing to DifferentiationInterface.jl PR, another bug from NonlinearSolve.jl has been spotted: Access to undefined reference when solving problems with Dual of BigFloat and fixed by generalizing the nonlinear solving with dual number to every package: refactor: Move dual nonlinear solving to NonlinearSolveBase

3. Implementing infinite final time boundary value problems solver, check out how BVPSUITE have managed to done it.

More TODOs:

1. Complete the PR for refactoring BoundaryValueDiffEq.jl from SparseDiffTools.jl to DifferentiationInterface.jl

2. Start working on the implementation of infinite final time BVP solver, and will start on this when the above feature is done.