Mesh adaptation error

[mmp14]

Hi,

I am running a Periodic Orbit continuation but it is failing with this error:

┌ Error: [Mesh-adaptation]. The mesh is non monotonic! Please report the error to the website of BifurcationKit.jl
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/RaJtn/src/periodicorbit/PeriodicOrbitCollocation.jl:1148

I cannot find where to report it on the website so I am posting here instead.

Here is the code I am using:

using Revise, Plots
using DifferentialEquations
using BifurcationKit

const BK = BifurcationKit

# # Define the function σ
function σ(v, v0, ϕ0 = 2.5, r=0.56)
    2 * ϕ0 / (1 + exp(r * (v0 - v)))
end

function lamnn(z, param)
    (;A_a, A_gs, A_gf, a_a, a_gs, a_gf, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, C12, C13, v0, v0_p2, ϕ0, r, ϕe1, ϕe2) = param
    y1, y6, y2, y7, y3, y8, y4, y9, y5, y10 = z
    [
        y6
        a_a * A_a * σ(C1 * y2 + C2 * y3 + C3 * (A_a / a_a) * ϕe1 + C11 * y4, v0) - 2 * a_a * y6 - a_a^2 * y1
        y7
        a_a * A_a * σ(C4 * y1, v0) - 2 * a_a * y7 - a_a^2 * y2
        y8
        a_gs * A_gs * σ(C5 * y1, v0) - 2 * a_gs * y8 - a_gs^2 * y3
        y9
        a_a * A_a * σ(C6 * y4 + C7 * y5 + C8 * (A_a / a_a) * ϕe2 + C12 * y1, v0_p2) - 2 * a_a * y9 - a_a^2 * y4
        y10
        a_gf * A_gf * σ(C9 * y4 + C10 * y5 + C13 * y1, v0) - 2 * a_gf * y10 - a_gf^2 * y5
    ]
end

z0 = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

par = (A_a = 3.25,
        A_gs = -22.0,
        A_gf = -30.0,
        a_a = 100.0,
        a_gs = 50.0,
        a_gf = 220.0, 
        C1 = 108.0,
        C2 = 33.7,
        C3 = 1.0,
        C4 = 135.0,
        C5 = 33.75,
        C6 = 70.0,
        C7 = 550.0,
        C8 = 1.0,
        C9 = 200.0,
        C10 = 100.0,
        C11 = 80.0,
        C12 = 200.0,
        C13 = 30.0,
        v0 = 6.0,
        v0_p2 = 1.0, 
        ϕ0 = 2.5,
        r = 0.56,
        ϕe1 = -100.0,
        ϕe2 = 0.0)

prob = BifurcationProblem(lamnn, z0, par, (@optic _.ϕe1);
	record_from_solution = (x, p; k...) -> (y1 = x[1]),)


optnewton = NewtonPar(tol = 1e-8, verbose = true)

 # # continuation options, we limit the parameter range for p
opts_br = ContinuationPar(p_min = -100.0, p_max = 500.0, dsmin = 0.000001, ds = 0.001, max_steps = 8000, newton_options = optnewton)

# # continuation of equilibria
br = continuation(prob, PALC(), opts_br;
	bothside = true, normC = norminf)

#Periodic orbit continuation
opts_br_po = ContinuationPar(p_min =-100.0, p_max = 400.0, max_steps = 4500, dsmin = 1e-7, ds = -1e-3, newton_options = NewtonPar(tol = 1e-8, verbose = true), detect_bifurcation = 3, tol_stability = 1e-3)

br_po = @time continuation(br, 4, opts_br_po,
                    PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
                    verbosity = 0,
                    alg = PALC(tangent = Bordered()),
                    linear_algo = COPBLS(),
                    normC = norminf,
                    callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
)

[rveltz]

Hi,

Thank you for posting this. I cannot reproduce. What is your version?

[mmp14]

I am using BifurcationKit v0.4.4

[mmp14]

I just upgraded to v0.4.6 and I still get the same error

[rveltz]

Can you please send the last continuation steps from the repl please?

[mmp14]

┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3198e-10 │ 0 │
└─────────────┴──────────────────────┴────────────────┘
┌ Warning: The precision on the Floquet multipliers is 0.15417195911976467.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease tol_stability in the option ContinuationPar or use a different method than FloquetColl for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/dbu6D/src/periodicorbit/Floquet.jl:535

┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3198e-10 │ 0 │
└─────────────┴──────────────────────┴────────────────┘
┌ Warning: The precision on the Floquet multipliers is 0.15417197535806146.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease tol_stability in the option ContinuationPar or use a different method than FloquetColl for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/dbu6D/src/periodicorbit/Floquet.jl:535

┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3198e-10 │ 0 │
└─────────────┴──────────────────────┴────────────────┘
┌ Warning: The precision on the Floquet multipliers is 0.1541719768370507.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease tol_stability in the option ContinuationPar or use a different method than FloquetColl for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/dbu6D/src/periodicorbit/Floquet.jl:535

┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3198e-10 │ 0 │
└─────────────┴──────────────────────┴────────────────┘
┌ Warning: The precision on the Floquet multipliers is 0.15417197036928565.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease tol_stability in the option ContinuationPar or use a different method than FloquetColl for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/dbu6D/src/periodicorbit/Floquet.jl:535

┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3198e-10 │ 0 │
└─────────────┴──────────────────────┴────────────────┘
┌ Warning: The precision on the Floquet multipliers is 0.15417196715802825.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease tol_stability in the option ContinuationPar or use a different method than FloquetColl for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/dbu6D/src/periodicorbit/Floquet.jl:535
┌ Error: [Mesh-adaptation]. The mesh is non monotonic! Please report the error to the website of BifurcationKit.jl

[mmp14]

┌─ Curve type: PeriodicOrbitCont from Hopf bifurcation point.
 ├─ Number of points: 3007
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter ϕe1 starts at 363.3126241517448, ends at 365.36578598365486
 ├─ Algo: PALC
 └─ Special points:

- #  1,       bp at ϕe1 ≈ +153.74586966 ∈ (+153.74572188, +153.74586966), |δp|=1e-04, [converged], δ = ( 1,  0), step = 1495
- #  2,       bp at ϕe1 ≈ +154.17680054 ∈ (+154.17680054, +154.45074524), |δp|=3e-01, [    guess], δ = (-1,  0), step = 1503
- #  3,       nd at ϕe1 ≈ +365.36578598 ∈ (+365.36578468, +365.36578598), |δp|=1e-06, [   guessL], δ = (-3, -2), step = 3006
- #  4, endpoint at ϕe1 ≈ +365.36578598,                                                                     step = 3006

[rveltz]

Pliease use ``` to put code

[rveltz]

Hi,

First. You can greatly improve the speed using

function lamnn!(dz, z, param, t = 0)
    (;A_a, A_gs, A_gf, a_a, a_gs, a_gf, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, C12, C13, v0, v0_p2, ϕ0, r, ϕe1, ϕe2) = param
    y1, y6, y2, y7, y3, y8, y4, y9, y5, y10 = z
    dz[1] = y6
    dz[2] = a_a * A_a * σ(C1 * y2 + C2 * y3 + C3 * (A_a / a_a) * ϕe1 + C11 * y4, v0) - 2 * a_a * y6 - a_a^2 * y1
    dz[3] = y7
    dz[4] = a_a * A_a * σ(C4 * y1, v0) - 2 * a_a * y7 - a_a^2 * y2
    dz[5] = y8
    dz[6] = a_gs * A_gs * σ(C5 * y1, v0) - 2 * a_gs * y8 - a_gs^2 * y3
    dz[7] = y9
    dz[8] = a_a * A_a * σ(C6 * y4 + C7 * y5 + C8 * (A_a / a_a) * ϕe2 + C12 * y1, v0_p2) - 2 * a_a * y9 - a_a^2 * y4
    dz[9] = y10
    dz[10] = a_gf * A_gf * σ(C9 * y4 + C10 * y5 + C13 * y1, v0) - 2 * a_gf * y10 - a_gf^2 * y5
    dz
end

z0 = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

prob = BifurcationProblem(lamnn!, z0, par, (@optic _.ϕe1);
    record_from_solution = (x, p; k...) -> (y1 = x[1]))

and then use:

PeriodicOrbitOCollProblem(100,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true)

[rveltz]

This seems to work fine for me

#Periodic orbit continuation
opts_br_po = ContinuationPar(p_min =-100.0, p_max = 1500.0, max_steps = 3000, dsmin = 0.000001, ds = -1e-2, dsmax = 0.1, newton_options = NewtonPar(tol = 1e-9, verbose = false), detect_bifurcation = 2, tol_stability = 1e-3, plot_every_step = 50)

args_po = (    record_from_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        _min, _max = extrema(xtt[1,:])
        return (max = _max,
                min = _min,
                amp = _max - _min,
                period = getperiod(p.prob, x, p.p))
    end,
    plot_solution = (x, p; k...) -> begin
        xtt = get_periodic_orbit(p.prob, x, p.p)
        arg = (marker = :d, markersize = 1)
        plot!(xtt.t, xtt[1:1,:]'; label = "y0", arg..., k...)
        plot!(br, subplot = 1)
        end,
    # we use the supremum norm
    normC = norminf)


br_po = @time continuation(br, 4, opts_br_po,
                    PeriodicOrbitOCollProblem(100,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true, K = 10, verbose_mesh_adapt = true, update_section_every_step = 10);
                    # ShootingProblem(15, prob_ode, Rodas5(); abstol = 1e-10, reltol = 1e-7, parallel = true);
                    verbosity = 1, plot = true,
                    δp = 0.2, 
                    alg = PALC(tangent = Bordered()),
                    linear_algo = COPBLS(),
                    normC = norminf,
                    callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
                    args_po...
            )

plot(br, br_po);
plot!(br_po, vars = (:param, :min))

[rveltz]

[rveltz]

Note that I incrase the discretization

[mmp14]

Thank you! I was able to get the same result as you. However, there are 2 fold bifurcations of limit cycles very close together around p = 160 that are not being detected. Is it possible to to that in BifurcationKit?

[rveltz]

Yes, they appear for PeriodicOrbitOCollProblem(60,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true, K = 10, verbose_mesh_adapt = true, update_section_every_step = 10)

[rveltz]

hence, for higher precision, they disappear. Where do you get that they should be here?

[mmp14]

This is from some previous analysis that was done in AUTO

[mmp14]

I see them now with PeriodicOrbitOCollProblem(60,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true, K = 10, verbose_mesh_adapt = true, update_section_every_step = 10) but they are marked as "bp" not "fold". Also the BP at p = 103 is a SNIC, but it also says "bp"

[rveltz]

A few things. You should use Auto the same discretization number.

In BK, we label the bifurcation points as bp in general. Going further requires the computation of normal forms. For folds however, there is a test function that can be used.

[mmp14]

Okay thanks. And what is a normal form and how do we compute them? What is the test function for folds?

[rveltz]

• You can pass event = BK.FoldDetectEvent, to continuation.
• For normal forms, see here and many others.
  AUTO detects SNIC?

[mmp14]

Re: AUTO, apparently it doesnt detect SNICs straight away, my bad!

I tried event = BK.FoldDetectEvent,, but the result seems the same

[rveltz]

can you share the auto code out of curiosity?

[rveltz]

opts_br_po = ContinuationPar(p_min =-100.0, p_max = 1500.0, max_steps = 3000, dsmin = 0.000001, ds = -1e-2, dsmax = 0.1, newton_options = NewtonPar(tol = 1e-9, verbose = false), detect_bifurcation = 1, tol_stability = 1e-3, plot_every_step = 100, detect_event = 2)

br_po = @time continuation(br, 4, opts_br_po,
                    PeriodicOrbitOCollProblem(60,4; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true, K = 100, verbose_mesh_adapt = false, update_section_every_step = 10);
                    verbosity = 3, plot = true,
                    δp = 0.2, 
                    alg = PALC(tangent = Bordered()),
                    linear_algo = COPBLS(),
                    event = BK.PairOfEvents(BK.FoldDetectEvent, BK.BifDetectEvent),
                    normC = norminf,
                    callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
                    args_po...
            )