I was very fortunate to be selected for this year’s SciML GSoC under NumFOCUS umbrella, a big thank you to my mentor Chris, Yingbo and Avik for letting me have this opportunity!!!

It has been a month so far, what have I been doing these days? Here is a summary.

Understand and compile MIRKDC program

After I finished my graduation project and thesis, I immediately started to understand the MIRKDC program, I learned the basic structure of MIRKDC and got the solver and driver programs running on my computer:

1 2	gfortran -o testing mirkdc.f mirkdc_driver.f -std=legacy blas_LINUX.a ./testing

Then we can see an interactive prompt where we can specify which method we want to use and the parameters of the given “Swirling flow III” problem.[^1] More details can be found in this blog post: https://erikqqy.github.io/2023/05/20/compile-mirkdc/

Develop defect control in BoundaryValueDiffEq.jl

Besides, I also tried to add a defect control routine to BoundaryValueDiffEq.jl. Right now, the basic structure is done, relating PR is here and the main flow can be summarized in a short description:

We compute the defect estimation in the whole interval, and generate new mesh to dimish the defect to get more accurate solution.

The detailed computation flow is depicted as follows:

Based on the discrete solution we already got using NLsolve.jl and MIRK tableau,

opt = isa(prob.problem_type, TwoPointBVProblem) ?
            alg.nlsolve(ConstructJacobian(jac_wrapper, vec_y), vec_y) :
            alg.nlsolve(ConstructJacobian(jac_wrapper, S, vec_y), vec_y)
nest_vector!(S.y, opt[1])

We can evaluate the weight polynomials(algorithm-dependent) at two sample points and to get the weights and their first derivatives by interp_weights. Then in every subinterval, we use interp_setup to compute the new stages in k_interp to prepare the construction of interpolants. We add the discrete solution, weights, weights derivative, discrete stages stored in k_discrete and additional stages stored in k_interp to evaluate the interpolant and compute the defect estimation in each subinterval. After we get the defect of the entire interval, we compute the defect norm simply by defect_norm = maximum(abs.(defect)). Then we estimate whether the defect_norm is bigger than the threshold(here, we set this threshold as 10%) and the tolerance.

After the estimation of the defect, all we need to do is to decide whether we need to refine the mesh by half_mesh(divide the original subinterval as two equal subintervals) or redistribute(generate a new mesh according to the defect on each subinterval to diminish the defect on each subinterval), this process takes place in mesh_selector. After the new mesh is generated, we restart the above process and reestimate the defect until the defect estimation is lower than the tolerance.

The detailed working flow can be visualized in the diagram below:

flow

Some things unclear

But there is still something I am unclear about. The most important one is the discrete stages of MIRK method in the entire interval, aka k_discrete. I know that in the MIRKDC routine, the discrete stages for MIRK method are computed in the process of computing PHI, which is the value of the residual function on each subinterval. This process has a correspondence part in BoundaryValueDiffEq.jl, that is the Φ! and update_K! part, where we update K in MIRK4GeneralCache during the nonlinear solving process.

1
2
3

struct MIRK4GeneralCache{kType} <: GeneralMIRKCache
    K::kType
end

function Φ!(S::BVPSystem{T}, TU::MIRKTableau, cache::AbstractMIRKCache) where {T}
    ####
    for i in 1:(N-1)
        update_K!(S, cache, TU, i, h)
    end
    ####
end

Where the update_K! is defined as:

function update_K!(S::BVPSystem, cache::AbstractMIRKCache, TU::MIRKTableau, i, h)
    M, N, residual, x, y, fun!, order = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order
    K, b = cache.K, TU.b
    c, v, X = TU.c, TU.v, TU.x

    function Kᵣ!(Kr, y, y₁, r)
        x_new = x[i] + c[r] * h
        y_new = (1 - v[r]) * y + v[r] * y₁
        if r > 1
            y_new += h * sum(j -> X[r, j] * K[j], 1:(r - 1))
        end
        fun!(Kr, y_new, S.p, x_new)
    end

    for r in 1:order
        Kᵣ!(K[r], y[i], y[i + 1], r)
    end
end

While I built the basic structure, there are still some things that haunted me for a long time. When I want to assemble all of the discrete stages K in all of the subintervals, the iteration seems weird, what I finally got is a Vector{Vector} which stored the same vector components. (And after the nonlinear solving process, the value of cache.K only stored the discrete stage in the last subinterval)

References

Ascher, Mattheij, and Russell, “Numerical Solution of Boundary Value Problems for Ordinary Differential Equations”, Classics in Applied Mathematics Series, Society for Industrial and Applied Mathematics, Philadelphia, 1995