Review of the past two weeks

In the last two weeks, my work can be summarized as two parts: Refactoring the error estimation, mesh refinement routines for COLDAE solvers and helping with the Lobatto IIIa-IIIc and RadauIIa methods for BVPs.

The previous work focus on the implementation and refactorization of discretized nonlinear solving stuff to get the accurate numerical result of the nonlinear collocation system, although the error estimation and mesh refinement part is implemented before, they still need further polishment to fix the potential bugs. With the recent pushes in collocation equations solving, I started on actually completing this part of work.

The basic workflow of error estimation and mesh refinement procedures can be summarized as flow:

# basic workflow

# basic model specifications
mesh1, f, g, J, dg = prob
sol1 = nonlinear_solving(mesh1)
adaptive && return

mesh2 = halve_mesh(mesh1)
if sol1.retcode == ReturnCode.Success
	sol2 = nonlinear_solving(mesh2) # solve on halved mesh
	err = error_estimate(sol1, sol2)
	if err > abstol
			mesh3 = mesh_selector(cache)
			info = ReturnCode.Failure # Start again
	end
else
	if length(mesh1) > max_num_subintervals
			info = ReturnCode.Failure
	else
			halve_mesh(mesh1)
			info = ReturnCode.Success # Restart
	end
end

It is similar with what we are doing in MIRK methods defect estimation and mesh refinement, except that we use the comparison of two solutions in different mesh to estimate the error norm.

Now the current implemented error estimate and mesh refinement routines are done, what’s left is the ForwardDiff stuff and sparse Jacobian configuation.

Helped with the completion of LobattoIIIa-LobattoIIIc and RadauIIa solvers. Finally all the expanded and nested version of FIRK solvers are done now, although the CI still have some issues to be addressed, it will be good to go in a few days.
Working on replacing solution objects when we are evaluating boundary conditions, after this PR is done, all the API in MIRK methods and Shooting methods would be unified(and they always should be), they would be changed to the interpolation in the solution in boundary conditions:

using BoundaryValueDiffEq
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -9.81 * sin(θ)
end
function bc!(residual, u, p, t)
    residual[1] = u(pi/4)[1] + pi / 2
    residual[2] = u(pi/2)[1] - pi / 2
end
prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
sol = solve(prob, MIRK4(), dt = 0.05)

This interpolation based indexing boundary conditions is better and we should change to this semantics long before😅.

More TODOs:

Refactor the implemented solvers to be compatible with the SciML style, and simplify the internal of Lobatto DAE solvers.
Address the current issues in FIRK solvers.