To make sure the project going on the right track with appropriate speed, I will record every week’s progress in this blog during the whole summer.

Review this week

It is really exciting that GSoC has been a week, in the first week, I opened a PR to add the defect control adaptivity for BoundaryValueDiffEq.jl, this PR provides the basic structure of the defect control technique for BVP solving, the detailed explanation can be seen in another blog post. Although it is still working in progress and there are still some problems in the computing process, I debugged the WIP PR and fixed some errors, I think the basic working flow is almost done. I noticed Avik had some work on the purge NLsolve.jl with NonlinearSolve.jl, wish this modification will make the solving process more accurate and make the discrete stages for MIRK algorithms easy to use.

While the previous blog post didn’t elaborate on the mesh refinement flow in detail, this week, I had a deeper look into the mesh refinement flow in MIRKDC, and here, I would like to record my understanding of how can we implement mesh refinement techniques in BoundaryValueDiffEq.jl and hopefully help me better understand the whole solving process.

The main mesh refinement part is located in the mesh_selector routine. First, we need to use the computed defect to generate mesh function values in the entire interval. The mesh function values of a subinterval can be computed using the below equation:

Then we can set r1 as the maximum of mesh function values, r2 as the weighted sum of all mesh function values, r3 as the average of weighted mesh function values. We compare the maximum of mesh function values r1(aka test if r1 is not sufficiently large, which means the distribution of mesh function values is sufficiently uniform), here we use r1 < rho*r3 to test r1 is not sufficiently large, if so, we halve the mesh. Otherwise(the mesh function values are sufficiently large, which means the distribution of mesh function values is sufficiently non-uniform), we redistribute a new mesh.

if r1 <= rho*r3
    ###
    mesh_new = half_mesh(mesh_current, n)
    ###
else
    ###
    mesh_new = redistribute(mesh_current, n, Nsub_star, s_hat)
    ###
end

Halve mesh

Halving mesh simply divides the subinterval equally as two equal subintervals.

function half_mesh(mesh_current::Vector, n::Int64)
    mesh_new = zeros(Float64, 2 * n + 1)
    mesh_new[1] = mesh_current[1]
    for i in 1:n
        mesh_new[2 * i + 1] = mesh_current[i + 1]
        mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2.0
    end
    return mesh_new
end

Redistribute mesh

The main idea of the mesh redistribution is to equidistribute the mesh function values.

We achieve the redistribution process by monitoring the total area under the mesh function Integral+Next_Piece, and comparing it with (the total area under the mesh function divided by Nsub_star(number of subintervals of the new mesh))

That is:

while k <= n
    next_piece = s_hat[k] * (mesh_current[k + 1] - t)
    if (integral + next_piece) > zeta
        mesh_new[i + 2] = (zeta - integral) / s_hat[k] + t
        t = mesh_new[i + 2]
        i = i + 1
        integral = 0
    else
        integral = integral + next_piece
        t = mesh_current[k + 1]
        k = k + 1
    end
end

Then we can get a new mesh with Nsub_star subintervals that satisfy our requirement.