[DifferentialEquations.jl] Solving Berman's equation

General Usage
1

Hi everyone,

I’m trying to solve the axisymetric Berman’s equation (see here) with DifferentialEquations.jl. This is a fourth order BVP :

\left[ \left( \frac{F}{r} \right) \left(\frac{F'}{r} \right)' - \left( \frac{F'}{r} \right)^2 \right]' - \frac{1}{R} \left[\frac{1}{r} \left( r \left( \frac{F'}{r}\right)' \right)' \right]' = 0,

with the associated boundary conditions:

\lim_{r \rightarrow 0} \left(\frac{F'}{r} \right)' = \lim_{r \rightarrow 0} \frac{F}{r} = F'(1) = 0,~ F(1) = 1.

R is a constant parameter corresponding to the flow Reynolds number. The solution is very close to F(r) = sin(\pi/2 r^2), especially when R \rightarrow +\infty.

I used to solve this equation using MATLAB bvp4c function. I am know trying DifferentialEquations.jl. I can recover the correct solution but as soon as I ask for a “fine” mesh, i.e \Delta r = 0.001, the computation time is very high (more than 10 minutes on my laptop) whereas it was 60s maximum for MATLAB bvp4c.

I think a clever initialization could help, but I didn’t find out how to specify the \sin(\pi/2r^2) solution as an initial solution.

Could anybody look at my implementation and give me some advice please? Especially:

  • what should look like the init. vector?
  • how to deal with the jacobian singularity in r=0?

Here is my code:

function f!(du,u,p,r)
    F0  = u[1]
    F1  = u[2]
    F2  = u[3]
    F3  = u[4]

    if (r < p.myeps)
        du .= [
            F1
            F2
            F3
            0.
        ]
    else
        du .= [
            F1
            F2
            F3
            1/r*( 2*F3-3*F2/r + 3*F1/(r*r) ) - p.Rey/r*( -F1*F2+F0*F3 -3*F0*F2/r + F1*F1/r +3*F0*F1/(r*r) )
        ]
    end
end

function bc!(residual, u, p, r)
    uleft = u[1]
    uright = u[end]
    residual .= [
        uright[1] - 1.
        uright[2]
        uleft[1]
        uleft[2]
    ]
end

rspan = (0.,1.)
p = (Rey = Rey, myeps = 1e-20)
u0 = [0., 0., 0., 0.] # I don't know how to use the init function
bvp = TwoPointBVProblem(f!, bc!, u0, rspan, p)
sol = solve(bvp, MIRK4(), dt=dt) # we need to use the MIRK4 solver for TwoPointBVProblem
2

I am not an expert on these particular kinds of problem, so maybe there is a bigger issue but I do see some unnecessary allocations in your derivative funciton for example:

Unless you have a particular reason not to do that, I would suggest you replace all these assignments with

        du .= F1,F2,F3,0.

which does not allocate anything. If I benchmark your f! function this way, I reduced the time by a factor of 2.

Taking this one step further, maybe it is worth to consider using StaticArrays since your state vector only has four elements. They are stack-allocated and are likely faster in your case as suggested by the DifferentialEquations.jl library:
https://diffeq.sciml.ai/stable/basics/faq/#My-ODE-is-solving-really-slow
https://diffeq.sciml.ai/stable/tutorials/ode_example/#ode_other_types

Ultimately, there may of course still be some issue that has to do with the BV problem in particular, I hope someone else who knows more about these chimes in in that case. I’d probably still at least try out the first option I’ve mentioned since its easily done.

3 Likes
3

Thank you. I’m bit ashamed that I did not present a code with SA, I was only focused on the package and not and my function. I will try that.

However I still believe that initializing properly the problem would greatly improve the computation time.

4

I have followed your advice. Here is the updated code:

function f!(du,u,p,r)
    F0  = u[1]
    F1  = u[2]
    F2  = u[3]
    F3  = u[4]

    if (r < p.myeps)
        du .= F1, F2, F3, 0.
    else
        du .= F1, F2, F3, 1/r*( 2*F3-3*F2/r + 3*F1/(r*r) ) - p.Rey/r*( -F1*F2+F0*F3 -3*F0*F2/r + F1*F1/r +3*F0*F1/(r*r) )
    end
end

function bc!(residual, u, p, r)
    uleft = u[1]
    uright = u[end]
    residual .= uright[1] - 1., uright[2], uleft[1], uleft[2]
end

rspan = (0.,1.)
p = (Rey = Rey, myeps = 1e-20)
u0 = [0., 0., 0., 0.] # I don't know how to use the init function
bvp = TwoPointBVProblem(f!, bc!, u0, rspan, p)
sol = solve(bvp, MIRK4(), dt=dt)

Using BenchmarkTools, I get the following results for the solve with dt=0.01 and Rey=2000.
Code from my first post: 2.291 s (39462640 allocations: 4.03 GiB)
Updated code : 2.229 s (35868878 allocations: 3.65 GiB)

So even if we’ve improved the performance of f, the solve performance is not significantly improved…

1 Like
5

Using ProfileView.jl, we can find out where all the runtime is coming from:

image
image1140×735 26 KB

All the runtime cost appears downstream of FiniteDiff.finite_diff_jacobian!, which is allocating tons of temporary arrays despite the availability of a cache for the Jacobian. @YingboMa might have some idea of what’s going wrong here.

4 Likes
6

Yeah, unfortunately, BoundaryValueDiffEq.jl doesn’t use optimized approximate Jacobian evaluations yet.

7

Have you tried a shooting method instead?

8

Is the shooting ok for TwoPointBVP? What algorithm do you recommend in association to the Shooting object?

And what should look like the init vector? The Shooting should be efficient only if I can provide an initial guess.

9

Looking at the sources, it seems that BVPSystem doesn’t handle user-provided jacobian function. Am I right?

10

Also, is not bv4pc adaptive?

11

Yeah, that’s right. Contributions are welcome. I wrote the code when I was in high school, so it’s very limited comparing to other software in SciML.