BifurcationKit: Mesh adaptation error

Specific Domains Modelling & Simulations
1

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
)
2

you can report here Issues · bifurcationkit/BifurcationKit.jl · GitHub

1 Like