search

bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
288 stars 34 forks source link

issue with user defined jacobian and codim2 continuation #152

Closed rveltz closed 1 month ago

rveltz commented 1 month ago

Hello,

I am using version “BifurcationKit v0.3.3” according to Julia. I agree my code does not make much sense, I am new to BifurcationKit.jl and was/am fumbling a bit.

I have removed the incorrect “jacobian_ma " and “jacobian” arguments that you pointed out. However, I now run into an additional problem – whenever I do not set " jacobian_ma” to anything, even when incorrect, I get the following error:

ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}}, Float64}, Float64, 5})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  Float64(::IrrationalConstants.Fourπ)
   @ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{…}, Float64}, Float64, 5})
    @ Base ./number.jl:7
  [2] setindex!
    @ ./array.jl:1024 [inlined]
  [3] JAC(u::SubArray{ForwardDiff.Dual{…}, 1, Vector{…}, Tuple{…}, true}, p::@NamedTuple{mu::ForwardDiff.Dual{…}, eta::Float64, k2::Float64})
    @ Main ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BailieKeane_2024_nonauto/Test:40
  [4] jacobian(pb::BifFunction{…}, x::SubArray{…}, p::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:72
  [5] jacobian(pb::BifurcationProblem{…}, x::SubArray{…}, p::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:234
  [6] (::FoldProblemMinimallyAugmented{…})(x::SubArray{…}, p::ForwardDiff.Dual{…}, params::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:35
  [7] (::FoldProblemMinimallyAugmented{…})(x::Vector{…}, params::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:48
  [8] (::BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}})(z::Vector{ForwardDiff.Dual{…}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172
  [9] vector_mode_dual_eval!(f::BifurcationKit.var"#352#353"{…}, cfg::ForwardDiff.JacobianConfig{…}, x::Vector{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
 [10] vector_mode_jacobian(f::BifurcationKit.var"#352#353"{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
 [11] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
 [12] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
 [13] jacobian
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172 [inlined]
 [14] _newton(prob::BifurcationKit.FoldMAProblem{…}, x0::Vector{…}, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(norminf), callback::typeof(BifurcationKit.cb_default), kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:88
 [15] _newton
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:62 [inlined]
 [16] newton
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:138 [inlined]
 [17] iterate(it::ContIterable{…}; _verbosity::Int64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:271
 [18] iterate
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250 [inlined]
 [19] continuation(it::ContIterable{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
 [20] continuation(prob::BifurcationKit.FoldMAProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::BorderingBLS{…}, bothside::Bool, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:539
 [21] continuation_fold(prob::BifurcationProblem{…}, alg::PALC{…}, foldpointguess::BorderedArray{…}, par::@NamedTuple{…}, lens1::Setfield.PropertyLens{…}, lens2::Setfield.PropertyLens{…}, eigenvec::Vector{…}, eigenvec_ad::Vector{…}, options_cont::ContinuationPar{…}; update_minaug_every_step::Int64, normC::typeof(norminf), bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, jacobian_ma::Symbol, compute_eigen_elements::Bool, usehessian::Bool, kind::BifurcationKit.FoldCont, record_from_solution::Nothing, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:497
 [22] continuation_fold
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:314 [inlined]
 [23] continuation_fold(prob::BifurcationProblem{…}, br::ContResult{…}, ind_fold::Int64, lens2::Setfield.PropertyLens{…}, options_cont::ContinuationPar{…}; alg::PALC{…}, normC::typeof(norminf), nev::Int64, start_with_eigen::Bool, bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:551
 [24] continuation_fold
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:511 [inlined]
 [25] #continuation#348
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/codim2.jl:237 [inlined]
 [26] top-level scope
    @ ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonautoBailieKeane_2024_nonauto/Test:195
 [27] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [28] top-level scope
    @ REPL[49]:1
in expression starting at /Users/johnbailie/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonautoBailieKeane_2024_nonauto/Test:195
Some type information was truncated. Use `show(err)` to see complete types.

On the other hand, when I set the “jacobian_ma = :a” for example, the correct branch of fold bifurcation is computed. I am honestly really lost on this, I believed I followed the documentation correctly?

Any assistance would be really appreciated.

FYI, this is my code

function JAC(u, p)
    m = p[1]
    eta = p[2]
    k2 = p[3]
    xn, yn, yl, xb = u

    ### Set main constants
    SV=10^6;                                                                               
    sig=2.1*10^4*SV;                                                                       
    alphaT=1.7*10^(-4);                                                                    
    alphaS=0.8*10^(-3);                                                                 
    ep=0.011;
    TaN=7;                                                                                 
    TaL=25                                                                              
    T0=2.65                                                                             
    S0 =35;                                                                                   
    SB =34.538;                                                                             
    yb = alphaS/alphaT * (SB - S0)/(TaN - T0);   
    g=3.17*10^-8;      
    VN=7.2106*10^15     #7*10^16;                                                                              
    VL=6.3515 * 10^16                                                                                
    VB=25 * VN     
    phi=alphaT*(TaN-T0);                                                                   
    xl=(TaL-T0)/(TaN-T0);    
    #                                                              
    TildeVL=VL/VN;                                                                         
    TildeVB=VB/VN;                                                                         
    K=5.456*SV/sig;    
    ps=phi*((yn-xn)-(yl-xl));   
    d = g*VN/sig

    # Functions                                                         
    k2n = k2;
    k2l = k2/3
    k1  = 10^-2*SV/sig;

    J = zeros(4,4)

    # FIRST ROW
    J[1,1] = (1/ep)*((0.5*k1*xb - 0.5*k2n*xb - 0.5*k1*xn + 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi^2*(-xb*xl - xl^2 + xb*xn + 3*xl*xn - 2*xn^2 + 
            xb*yl + xl*yl - 2*xn*yl - xb*yn - xl*yn + 2*xn*yn)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-d - K - 0.5*k1 - 
            0.5*k2n + (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(-xb - 3*xl + 4*xn + 2*yl - 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))

    J[1,2] = (1/ep)*((-0.5*k1*xb + 0.5*k2n*xb + 0.5*k1*xn - 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(phi*(xl^2 
            + xn*(2*xn + 2*yl - 2*yn) + xl*(-3*xn - yl + yn) + xb*(xl - xn - yl + yn)) + ep*(0.5*xb + 0.5*xl - xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))

    J[1,3] = -((phi*(xb + xl - 2*xn)*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(xl - xn - yl + yn) + ep*sinh((2*phi*(xl - xn - yl + yn))/ep)))/(2*ep))

    J[1,4] = 0.5*k1 + 0.5*k2n + (0.5*k1 - 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(xl - xn - yl + yn)*tanh((phi*(xl - xn - yl + yn))/ep)

    # SECOND ROW
    J[2,1] = (1/ep)*((0.5*k1*yb - 0.5*k2n*yb - 0.5*k1*yn + 0.5*k2n*yn)*sech((eta + xn - yn)/ep)^2 + phi*sech((phi*(xl - xn - yl 
            + yn))/ep)^2*(phi*(-xl*yb + xn*yb - xl*yl + xn*yl + yb*yl + yl^2 + 2*xl*yn - 2*xn*yn - yb*yn - 3*yl*yn + 2*yn^2) + ep*(-0.5*yb 
            - 0.5*yl + yn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))

    J[2,2] = (1/ep)*((-0.5*k1*yb + 0.5*k2n*yb + 0.5*k1*yn - 0.5*k2n*yn)*sech((eta + xn - yn)/ep)^2 + 
            phi^2*(xl*yb - xn*yb + xl*yl - xn*yl - yb*yl - yl^2 - 2*xl*yn + 2*xn*yn + yb*yn + 3*yl*yn - 
            2*yn^2)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-K - 0.5*k1 - 0.5*k2n + (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) + 
            phi*(-2*xl + 2*xn + yb + 3*yl - 4*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))

    J[2,3] = K + (1/(2*ep))*phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(yb + yl - 2*yn)*(-xl + xn + yl - yn) + 
            ep*(-yb + yn)*sinh((2*phi*(xl - xn - yl + yn))/ep)) + phi*(xl - xn - 2*yl + 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep)

    J[2,4] = 0.0

    # THIRD ROW
    J[3,1] = -((phi*(yb - 2*yl + yn)*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(xl - xn - yl + yn) + ep*sinh((2*phi*(xl - xn - yl + yn))/ep)))/(2*ep*TildeVL))

    J[3,2] = (1/TildeVL)*(K + (1/(2*ep))*phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(yb - 2*yl + yn)*(xl - xn - yl + yn) 
            + ep*(yb - yl)*sinh((2*phi*(xl - xn - yl + yn))/ep)) + phi*(xl - xn - 2*yl + 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep))

    J[3,3] = (1/(ep*TildeVL))*((-0.5*k1*yb + 0.5*k2l*yb + 0.5*k1*yl - 0.5*k2l*yl)*sech((eta + xl - yl)/ep)^2 + phi^2*(-xl*yb + xn*yb + 2*xl*yl - 2*xn*yl + yb*yl - 
            2*yl^2 - xl*yn + xn*yn - yb*yn + 3*yl*yn - yn^2)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-K - 0.5*k1 - 0.5*k2l + (-0.5*k1 
            + 0.5*k2l)*tanh((eta + xl - yl)/ep) + phi*(-2*xl + 2*xn - yb + 4*yl - 3*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))

    J[3,4] = 0.0

    # FOURTH ROW
    J[4,1] = (1/(ep*TildeVB))*((-0.5*k1*xb + 0.5*k2n*xb + 0.5*k1*xn - 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi^2*(-xl^2 + xl*(yl - yn) 
            + xn*(xn + yl - yn) + xb*(2*xl - 2*xn - 2*yl + 2*yn))*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(0.5*k1 + 0.5*k2n + (0.5*k1 
            - 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(2*xb - 2*xn - yl + yn)*tanh((phi*(xl - xn - yl + yn))/ep)))

    J[4,2] = (1/(ep*TildeVB))*((0.5*k1*xb - 0.5*k2n*xb - 0.5*k1*xn + 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 
            + phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(phi*(xl^2 - xn^2 - xl*yl - xn*yl + xb*(-2*xl + 2*xn 
            + 2*yl - 2*yn) + xl*yn + xn*yn) + ep*(-xb + 0.5*xl + 0.5*xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))

    J[4,3] = (1/(ep*TildeVB))*((0.5*k1*xb - 0.5*k2l*xb - 0.5*k1*xl + 0.5*k2l*xl)*sech((eta + xl - yl)/ep)^2 + phi*sech((phi*(xl - xn - yl 
            + yn))/ep)^2*(phi*(-xl^2 + xn^2 + xl*yl + xn*yl - xl*yn - xn*yn + xb*(2*xl - 2*xn - 2*yl + 2*yn)) + ep*(xb 
            - 0.5*xl - 0.5*xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))

    J[4,4] = (1/TildeVB)*(-k1 - 0.5*k2l - 0.5*k2n + (-0.5*k1 + 0.5*k2l)*tanh((eta + xl - yl)/ep) 
            + (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(-2*xl + 2*xn + 2*yl - 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep))

    # println(J)
    return J 
end

function ODEHandle2!(du, u, p)
    m = p[1]
    eta = p[2]
    k2 = p[3]
    xn, yn, yl, xb = u
    ### Set main constants
    SV=10^6;                                                                               
    sig=2.1*10^4*SV;                                                                       
    alphaT=1.7*10^(-4);                                                                    
    alphaS=0.8*10^(-3);                                                                 
    ep=0.011;
    TaN=7;                                                                                 
    TaL=25                                                                              
    T0=2.65                                                                             
    S0 =35;                                                                                   
    SB =34.538;                                                                             
    yb = alphaS/alphaT * (SB - S0)/(TaN - T0);   
    g=3.17*10^-8;      
    VN=7.2106*10^15     #7*10^16;                                                                              
    VL=6.3515 * 10^16                                                                                
    VB=25 * VN     
    phi=alphaT*(TaN-T0);                                                                   
    xl=(TaL-T0)/(TaN-T0);    
    #                                                              
    TildeVL=VL/VN;                                                                         
    TildeVB=VB/VN;                                                                         
    K=5.456*SV/sig;    
    ps=phi*((yn-xn)-(yl-xl));   
    d = g*VN/sig

    # Functions                                                         
    k2n = k2;
    k2l = k2/3
    k1  = 10^-2*SV/sig;

    kappaN=k1+0.5*(k2n-k1)*(1+tanh((yn-xn-eta)/ep));                                        
    kappaL=k1+0.5*(k2l-k1)*(1+tanh((yl-xl-eta)/ep));                                                                                                                                                      
    PsiPlus=1/2*(1+tanh(ps/ep))*ps                ;                                             
    PsiMinus=1/2*(1-tanh(ps/ep))*ps               ;                                            
    Psi=ps*tanh(ps/ep)                            ;                            
    # Set up function output                                                                  
    du[1] = -d*(xn-1)+(K+Psi)*(xl-xn)+(PsiMinus-PsiPlus-kappaN)*(xn-xb)    
    du[2] = m+(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaN)*(yn-yb)          
    du[3] = 1/TildeVL*(-m-(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaL)*(yl-yb))
    du[4] = 1/TildeVB*((kappaN-PsiMinus+PsiPlus)*(xn-xb)+(kappaL-PsiMinus+PsiPlus)*(xl-xb))

    du 
end

# options
opts = ContinuationPar(
        p_min = -3.0,
        p_max = 0.01,
        detect_bifurcation = 3,
        nev = 4,
        max_steps = 50000,
        dsmin = 1e-15,
        ds = 1e-3,
        dsmax = 2e-3,
        n_inversion = 8,
        max_bisection_steps =150,
        detect_event = 0
)

# Equilibria
z0 = [1.242167, -1.260954, 0.2613593, 3.190049]

# Define the bifurcation problem
prob = BifurcationProblem(
    ODEHandle2!,
    z0,
    par_tm,
    (@lens _.mu),
    J=JAC)

# Codimension-one continuation 
br = continuation(
    prob,
    PALC(),
    J=JAC,
    opts;
    normC = norminf,
    bothside=true,
    detect_bifurcation=3)

# Codimension-two continuation from the branch point at index 2
indFold = 2
sn_codim2 = continuation(br, indFold, (@lens _.eta), J=JAC,
    opts;
    bothside=true,
    normC = norminf,
    bdlinsolver = MatrixBLS(),
    update_minaug_every_step = 1,
    detect_codim2_bifurcation = 2,
    detect_loop=1
)
rveltz commented 1 month ago

The code is missing par_tm, can you provide it please?

ShurimasKitten commented 1 month ago

it should be this,

par_tm = (mu = -37.5*SV/sig, eta = -3.0, k2 = 8.0*SV/sig)

rveltz commented 1 month ago

But I dont have SV and sig to build par_tm

rveltz commented 1 month ago

I think your problem comes from your JAC function because your J is a Matrix{Float64}. Try not passing J as an argument to any function.

rveltz commented 1 month ago 

function ODEHandle2!(du, u, p)
    (;mu, eta, k2) = p
    m = mu
    xn, yn, yl, xb = u
    ### Set main constants
    SV=10^6;                                                                               
    sig=2.1*10^4*SV;                                                                       
    alphaT=1.7*10^(-4);                                                                    
    alphaS=0.8*10^(-3);                                                                 
    ep=0.011;
    TaN=7;                                                                                 
    TaL=25                                                                              
    T0=2.65                                                                             
    S0 =35;                                                                                   
    SB =34.538;                                                                             
    yb = alphaS/alphaT * (SB - S0)/(TaN - T0);   
    g=3.17*10^-8;      
    VN=7.2106*10^15     #7*10^16;                                                                              
    VL=6.3515 * 10^16                                                                                
    VB=25 * VN     
    phi=alphaT*(TaN-T0);                                                                   
    xl=(TaL-T0)/(TaN-T0);    
    #                                                              
    TildeVL=VL/VN;                                                                         
    TildeVB=VB/VN;                                                                         
    K=5.456*SV/sig;    
    ps=phi*((yn-xn)-(yl-xl));   
    d = g*VN/sig

    # Functions                                                         
    k2n = k2;
    k2l = k2/3
    k1  = 10^-2*SV/sig;

    kappaN=k1+0.5*(k2n-k1)*(1+tanh((yn-xn-eta)/ep));                                        
    kappaL=k1+0.5*(k2l-k1)*(1+tanh((yl-xl-eta)/ep));                                                                                                                                                      
    PsiPlus=1/2*(1+tanh(ps/ep))*ps                ;                                             
    PsiMinus=1/2*(1-tanh(ps/ep))*ps               ;                                            
    Psi=ps*tanh(ps/ep)                            ;                            
    # Set up function output                                                                  
    du[1] = -d*(xn-1)+(K+Psi)*(xl-xn)+(PsiMinus-PsiPlus-kappaN)*(xn-xb)    
    du[2] = m+(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaN)*(yn-yb)          
    du[3] = 1/TildeVL*(-m-(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaL)*(yl-yb))
    du[4] = 1/TildeVB*((kappaN-PsiMinus+PsiPlus)*(xn-xb)+(kappaL-PsiMinus+PsiPlus)*(xl-xb))

    du 
end

# options
opts = ContinuationPar(
        p_min = -3.0,
        p_max = 0.01,
        detect_bifurcation = 3,
        nev = 4,
        max_steps = 50000,
        dsmin = 1e-15,
        ds = 1e-3,
        dsmax = 2e-3,
        n_inversion = 4,
        max_bisection_steps = 150,
        detect_event = 0
)

# Equilibria
z0 = [1.242167, -1.260954, 0.2613593, 3.190049]

SV=10^6;                                                                               
sig=2.1*10^4*SV; 
par_tm = (mu = -37.5*SV/sig, eta = -3.0, k2 = 8.0*SV/sig)

# Define the bifurcation problem
prob = BifurcationProblem(
    ODEHandle2!,
    z0,
    par_tm,
    (@lens _.mu),
#     J=JAC
    )

# Codimension-one continuation 
br = continuation(
    prob,
    PALC(),
#     J=JAC,
    opts;
    normC = norminf,
    bothside=true,
    # this is bad! detect_bifurcation is only passed to ContinuationPar
#     detect_bifurcation=3
    )
ShurimasKitten commented 1 month ago

But I dont have SV and sig to build par_tm

I am so sorry! That was a clueless oversight by me.

SV = 10^6 sig = 2.1 10^4SV par_tm = (mu = -37.5SV/sig, eta = -3.0, k2 = 8.0SV/sig)

Thank you so much for your reply though! I really appreciate it and your comments :) I will try what you suggests

rveltz commented 1 month ago

Another thing, you have very different parameter values with huge ranges. This will imped numerical continuation but also your simulations, etc. Try using adimensionalized values.

rveltz commented 1 month ago

you should also do this 

SV=10^6;                                                                               
sig=2.1*10^4*SV; 
par_tm = (mu = -37.5*SV/sig, eta = -3.0, k2 = 8.0*SV/sig, SV = SV, sig = sig)

function ODEHandle2!(du, u, p)
    (;mu, eta, k2, SV, sig) = p
    m = mu
    xn, yn, yl, xb = u
    ### Set main constants

    alphaT=1.7*10^(-4);                                                                    
    alphaS=0.8*10^(-3);                                                                 
    ep=0.011;
    TaN=7;                                                                                 
    TaL=25                                                                              
    T0=2.65                                                                             
    S0 =35;                                                                                   
    SB =34.538;                                                                             
    yb = alphaS/alphaT * (SB - S0)/(TaN - T0);   
    g=3.17*10^-8;      
    VN=7.2106*10^15     #7*10^16;                                                                              
    VL=6.3515 * 10^16                                                                                
    VB=25 * VN     
    phi=alphaT*(TaN-T0);                                                                   
    xl=(TaL-T0)/(TaN-T0);    
    #                                                              
    TildeVL=VL/VN;                                                                         
    TildeVB=VB/VN;                                                                         
    K=5.456*SV/sig;    
    ps=phi*((yn-xn)-(yl-xl));   
    d = g*VN/sig

    # Functions                                                         
    k2n = k2;
    k2l = k2/3
    k1  = 10^-2*SV/sig;

    kappaN=k1+0.5*(k2n-k1)*(1+tanh((yn-xn-eta)/ep));                                        
    kappaL=k1+0.5*(k2l-k1)*(1+tanh((yl-xl-eta)/ep));                                                                                                                                                      
    PsiPlus=1/2*(1+tanh(ps/ep))*ps                ;                                             
    PsiMinus=1/2*(1-tanh(ps/ep))*ps               ;                                            
    Psi=ps*tanh(ps/ep)                            ;                            
    # Set up function output                                                                  
    du[1] = -d*(xn-1)+(K+Psi)*(xl-xn)+(PsiMinus-PsiPlus-kappaN)*(xn-xb)    
    du[2] = m+(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaN)*(yn-yb)          
    du[3] = 1/TildeVL*(-m-(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaL)*(yl-yb))
    du[4] = 1/TildeVB*((kappaN-PsiMinus+PsiPlus)*(xn-xb)+(kappaL-PsiMinus+PsiPlus)*(xl-xb))

    du 
end
ShurimasKitten commented 1 month ago

Hello,

I think you are correct that J=JAC might be the problem here. Running the program with 

prob = BifurcationProblem(
    NonAutoAdvec!,
    #J = JAC, 
    z0,
    par_tm,
    (@lens _.mu),
    )

br = continuation(
    prob,
    PALC(),
    opts;
    normC = norminf,
    bothside=true,
    )

indFold = 3
sn_codim2 = continuation(br, indFold, (@lens _.eta),
    opts;
    bothside=true,
    normC = norminf,
    bdlinsolver = MatrixBLS(),
    update_minaug_every_step = 1,
)

results in the correct the bifurcation diagram.

However, when using the follow bifurcation problem 

prob = BifurcationProblem(
    ODEHandle2!,
    J = JAC, 
    z0,
    par_tm,
    (@lens _.mu),
    )

one gets the following error + stack trace 

julia> include("Test")
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{…}, Float64}, Float64, 5})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  Float64(::IrrationalConstants.Fourπ)
   @ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{…}, Float64}, Float64, 5})
    @ Base ./number.jl:7
  [2] setindex!
    @ ./array.jl:1024 [inlined]
  [3] JAC(u::SubArray{…}, p::@NamedTuple{…})
    @ Main ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BailieKeane_2024_nonauto/JAC.jl:38
  [4] jacobian(pb::BifFunction{…}, x::SubArray{…}, p::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:72
  [5] jacobian(pb::BifurcationProblem{…}, x::SubArray{…}, p::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:234
  [6] (::FoldProblemMinimallyAugmented{…})(x::SubArray{…}, p::ForwardDiff.Dual{…}, params::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:35
  [7] (::FoldProblemMinimallyAugmented{…})(x::Vector{…}, params::@NamedTuple{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:48
  [8] (::BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}})(z::Vector{ForwardDiff.Dual{…}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172
  [9] vector_mode_dual_eval!(f::BifurcationKit.var"#352#353"{…}, cfg::ForwardDiff.JacobianConfig{…}, x::Vector{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
 [10] vector_mode_jacobian(f::BifurcationKit.var"#352#353"{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
 [11] jacobian(f::Function, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…}, ::Val{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
 [12] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
 [13] jacobian
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172 [inlined]
 [14] _newton(prob::BifurcationKit.FoldMAProblem{…}, x0::Vector{…}, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(norminf), callback::typeof(BifurcationKit.cb_default), kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:88
 [15] _newton
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:62 [inlined]
 [16] newton
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:138 [inlined]
 [17] iterate(it::ContIterable{…}; _verbosity::Int64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:271
 [18] iterate
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250 [inlined]
 [19] continuation(it::ContIterable{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
 [20] continuation(prob::BifurcationKit.FoldMAProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::BorderingBLS{…}, bothside::Bool, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:539
 [21] continuation_fold(prob::BifurcationProblem{…}, alg::PALC{…}, foldpointguess::BorderedArray{…}, par::@NamedTuple{…}, lens1::Setfield.PropertyLens{…}, lens2::Setfield.PropertyLens{…}, eigenvec::Vector{…}, eigenvec_ad::Vector{…}, options_cont::ContinuationPar{…}; update_minaug_every_step::Int64, normC::typeof(norminf), bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, jacobian_ma::Symbol, compute_eigen_elements::Bool, usehessian::Bool, kind::BifurcationKit.FoldCont, record_from_solution::Nothing, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:497
 [22] continuation_fold
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:314 [inlined]
 [23] continuation_fold(prob::BifurcationProblem{…}, br::ContResult{…}, ind_fold::Int64, lens2::Setfield.PropertyLens{…}, options_cont::ContinuationPar{…}; alg::PALC{…}, normC::typeof(norminf), nev::Int64, start_with_eigen::Bool, bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, kwargs::@Kwargs{…})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:551
 [24] continuation_fold
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:511 [inlined]
 [25] #continuation#348
    @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/codim2.jl:237 [inlined]
 [26] top-level scope
    @ ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonauto/Test:44
 [27] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [28] top-level scope
    @ REPL[21]:1
in expression starting at /Users/johnbailie/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonauto/Test:44
Some type information was truncated. Use `show(err)` to see complete types.

So it seems that the Jacobian is indeed the issue here. Do you have an idea for a fix for this? Would changing J from a Matrix{Float64} to another data type be a starting point? It would be really nice to use the Jacobian during the continuation ^.^

Also, on your remark about non-dimensionalisation. The model is already non-dimensional, and the non-dimensional parameters in the model i.e. d, TildeVL, TildeVB, K, psi etc end up being quite reasonable values even though some of their dimensional dependences are a bit insane :) Climate models are a bit funky like that ^.^

rveltz commented 1 month ago

Yes I know! It is the line :

J = zeros(4,4)

is the problem. Change it to:

  J = similar(u, 4,4)

For automatic differentiation, the type of x changes to ForwardDiff.Dual. And when you write

J[1,1] = x[1]

you are trying to fit a ForwardDiff.Dual into Float64

rveltz commented 1 month ago

I am changing the title because the error stems from another reason

ShurimasKitten commented 1 month ago

With this I was able to solve the original problem I wanted to solve and the bifurcation diagram is drawn as expected.

Thank you so much!!! I would never have spotted that 🫠

rveltz commented 1 month ago

Should we close this?

The original problem was about not detecting user passed events in codim 2 continuation

ShurimasKitten commented 1 month ago

You can close it :)