Skip to content

Commit

Permalink
remove Setfield from tests
Browse files Browse the repository at this point in the history
  • Loading branch information
rveltz committed Sep 15, 2024
1 parent 6ac3869 commit 5b957b2
Show file tree
Hide file tree
Showing 20 changed files with 92 additions and 83 deletions.
20 changes: 12 additions & 8 deletions src/periodicorbit/FlowDE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,11 @@ function Flow(prob::Union{ODEProblem, EnsembleProblem, DAEProblem}, alg; kwargs.
return FlowDE(prob, alg, nothing, nothing, kwargs, get(kwargs, :callback, nothing), nothing, nothing, 1e-8)
end

function Flow(prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2; kwargs...)
function Flow(prob1::Union{ODEProblem, EnsembleProblem},
alg1,
prob2::Union{ODEProblem, EnsembleProblem},
alg2;
kwargs...)
return FlowDE(prob1, alg1, prob2, alg2, kwargs, get(kwargs, :callback, nothing), nothing, nothing, 1e-8)
end
####################################################################################################
Expand All @@ -66,7 +70,7 @@ end
function _flow(x, p, tm, pb::ODEProblem, alg; kwargs...)
_prob = remake(pb; u0 = x, tspan = (zero(tm), tm), p = p)
# the use of concrete_solve makes it compatible with Zygote
sol = solve(_prob, alg; save_everystep = false, kwargs...)
sol = SciMLBase.solve(_prob, alg; save_everystep = false, kwargs...)
return (t = sol.t[end], u = sol.u[end])
end
####################################################################################################
Expand All @@ -82,7 +86,7 @@ function evolve(fl::FlowDE{T1}, x::AbstractArray, p, tm; kw...) where {T1 <: Ens
# see docs at https://docs.sciml.ai/dev/features/ensemble/#Performing-an-Ensemble-Simulation-1
_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(eltype(tm[ii])), tm[ii]), p = p)
_epb = setproperties(fl.prob, output_func = (sol, i) -> ((t = sol.t[end], u = sol.u[end]), false), prob_func = _prob_func)
sol = solve(_epb, fl.alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, fl.kwargsDE..., kw...)
sol = SciMLBase.solve(_epb, fl.alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, fl.kwargsDE..., kw...)
# sol.u contains a vector of tuples (sol_i.t[end], sol_i[end])
return sol.u
end
Expand All @@ -92,13 +96,13 @@ function dflowMonoSerial(x::AbstractVector, p, dx, tm, pb::ODEProblem, alg; k...
n = length(x)
_prob = remake(pb; u0 = vcat(x, dx), tspan = (zero(tm), tm), p = p)
# the use of concrete_solve makes it compatible with Zygote
sol = solve(_prob, alg, save_everystep = false; k...)[end]
sol = SciMLBase.solve(_prob, alg, save_everystep = false; k...)[end]
return (t = tm, u = sol[1:n], du = sol[n+1:end])
end

function dflow_fdSerial(x, p, dx, tm, pb::ODEProblem, alg; δ = convert(eltype(x), 1e-9), kwargs...)
sol1 = _flow(x .+ δ .* dx, p, tm, pb, alg; kwargs...).u
sol2 = _flow(x , p, tm, pb, alg; kwargs...).u
sol2 = _flow(x , p, tm, pb, alg; kwargs...).u
return (t = tm, u = sol2, du = (sol1 .- sol2) ./ δ)
end

Expand All @@ -118,7 +122,7 @@ function jvp(fl::FlowDE{T1}, x::AbstractArray, p, dx, tm; kw...) where {T1 <: E
N = size(x, 1)
_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = vcat(x[:, ii], dx[:, ii]), tspan = (zero(tm[ii]), tm[ii]), p = p)
_epb = setproperties(fl.probMono, output_func = (sol,i) -> ((t = sol.t[end], u = sol[end][1:N], du = sol[end][N+1:end]), false), prob_func = _prob_func)
sol = solve(_epb, fl.algMono, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kw...)
sol = SciMLBase.solve(_epb, fl.algMono, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kw...)
return sol.u
end

Expand All @@ -137,13 +141,13 @@ end
# this function takes into account a parameter passed to the vector field and returns the full solution from the ODE solver. This is useful in Poincare Shooting to extract the period.
function evolve(fl::FlowDE{T1}, ::Val{:Full}, x::AbstractArray, p, tm; kw...) where {T1 <: ODEProblem}
_prob = remake(fl.prob; u0 = x, tspan = (zero(tm), tm), p = p)
sol = solve(_prob, fl.alg; fl.kwargsDE..., kw...)
sol = SciMLBase.solve(_prob, fl.alg; fl.kwargsDE..., kw...)
end

function evolve(fl::FlowDE{T1}, ::Val{:Full}, x::AbstractArray, p, tm; kw...) where {T1 <: EnsembleProblem}
_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(eltype(tm[ii])), tm[ii]), p = p)
_epb = setproperties(fl.prob, prob_func = _prob_func)
sol = solve(_epb, fl.alg, EnsembleThreads(); trajectories = size(x, 2), fl.kwargsDE..., kw...)
sol = SciMLBase.solve(_epb, fl.alg, EnsembleThreads(); trajectories = size(x, 2), fl.kwargsDE..., kw...)
end

function evolve(fl::FlowDE{T1}, ::Val{:SerialTimeSol}, x::AbstractArray, par, δt; k...) where {T1 <: ODEProblem}
Expand Down
21 changes: 13 additions & 8 deletions src/periodicorbit/ShootingDE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,12 @@ end
ShootingProblem(M::Int, prob::ODEType, alg; kwargs...) = ShootingProblem(prob, alg, M, nothing; kwargs...)

# idem but with an ODEProblem to define the derivative of the flow
function ShootingProblem(prob1::ODEType, alg1, prob2::ODEType, alg2, ds, section; parallel = false, par = prob1.p, kwargs...)
function ShootingProblem(prob1::ODEType, alg1,
prob2::ODEType, alg2,
ds, section;
parallel = false,
par = prob1.p,
kwargs...)
_M = length(ds)
parallel = _M == 1 ? false : parallel
_pb1 = parallel ? EnsembleProblem(prob1) : prob1
Expand All @@ -67,13 +72,13 @@ end
### POINCARE SHOOTING
####################################################################################################
function PoincareShootingProblem(prob::ODEProblem,
alg,
hyp::SectionPS;
δ = 1e-8,
interp_points = 50,
parallel = false,
par = prob.p,
kwargs...)
alg,
hyp::SectionPS;
δ = 1e-8,
interp_points = 50,
parallel = false,
par = prob.p,
kwargs...)
pSection(out, u, t, integrator) = (hyp(out, u); out .*= integrator.iter > 1)
affect!(integrator, idx) = terminate!(integrator)
# we put nothing option to have an upcrossing
Expand Down
2 changes: 1 addition & 1 deletion test/COModel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ br = continuation(prob, PALC(), opts_br; normC = norminf, bothside = true)
@test br.specialpoint[4].param 1.04204851
@test br.specialpoint[5].param 1.05158367
####################################################################################################
@set! opts_br.newton_options = NewtonPar(max_iterations = 10, tol = 1e-12)
@reset opts_br.newton_options = NewtonPar(max_iterations = 10, tol = 1e-12)

sn_codim2 = continuation(br, 3, (@optic _.k),
ContinuationPar(opts_br, p_max = 2.2, p_min = 0., ds = -0.001, dsmax = 0.05, n_inversion = 8, max_steps = 50) ;
Expand Down
10 changes: 5 additions & 5 deletions test/codim2.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# using Revise
using BifurcationKit
using Test, Setfield, LinearAlgebra
using Test, LinearAlgebra
# using Plots
const BK = BifurcationKit
####################################################################################################
Expand All @@ -22,7 +22,7 @@ prob = BifurcationProblem(COm, z0, par_com, (@optic _.q2); record_from_solution

opts_br = ContinuationPar(p_min = 0.6, p_max = 2.5, ds = 0.002, dsmax = 0.01, n_inversion = 4, detect_bifurcation = 3, max_bisection_steps = 25, nev = 2, max_steps = 20000)

# @set! opts_br.newton_options.verbose = true
# @reset opts_br.newton_options.verbose = true
alg = PALC()
br = @time continuation(prob, alg, opts_br;
plot = false, verbosity = 0, normC = norminf,
Expand Down Expand Up @@ -51,8 +51,8 @@ snpt = get_normal_form(br, 3)
@test snpt.nf.b2 0.2693795490512864 rtol = 1e-3
@test snpt.nf.b3 12.340786210650833 rtol = 1e-3

@set! opts_br.newton_options.verbose = false
@set! opts_br.newton_options.max_iterations = 10
@reset opts_br.newton_options.verbose = false
@reset opts_br.newton_options.max_iterations = 10

sn = newton(br, 3; options = opts_br.newton_options, bdlinsolver = MatrixBLS())
# printstyled(color=:red, "--> guess for SN, p = ", br.specialpoint[2].param, ", psn = ", sn[1].p)
Expand Down Expand Up @@ -96,7 +96,7 @@ hppt = get_normal_form(br, 2)
@test hppt.nf.a 2.546719962189168 + 1.6474887797814664im
@test hppt.nf.b 4.3536804635557855 + 15.441272421860365im

@set! opts_br.newton_options.verbose = false
@reset opts_br.newton_options.verbose = false

hp = BK.newton_hopf(br, 2; options = opts_br.newton_options, start_with_eigen = true)
# printstyled(color=:red, "--> guess for HP, p = ", br.specialpoint[1].param, ", php = ", hp[1].p)
Expand Down
10 changes: 5 additions & 5 deletions test/codim2PO-OColl.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,9 @@ using OrdinaryDiffEq
prob_de = ODEProblem(Pop!, z0, (0,600.), par_pop)
alg = Rodas5()
# alg = Vern9()
sol = solve(prob_de, alg)
sol = OrdinaryDiffEq.solve(prob_de, alg)
prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-8, abstol = 1e-10)
sol = solve(prob_de, Rodas5())
sol = OrdinaryDiffEq.solve(prob_de, Rodas5())
################################################################################
argspo = (record_from_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, p.p)
Expand Down Expand Up @@ -95,8 +95,8 @@ pd_po_coll = continuation(brpo_pd, 1, (@optic _.b0), opts_pocoll_pd;
################################################################################
# find the NS case
par_pop2 = @set par_pop.b0 = 0.4
sol2 = solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = solve(remake(sol2.prob, tspan = (0, 10), u0 = sol2[end]), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(sol2.prob, tspan = (0, 10), u0 = sol2[end]), Rodas5())

probcoll, ci = generate_ci_problem(PeriodicOrbitOCollProblem(26, 3; update_section_every_step = 0), re_make(prob, params = sol2.prob.p), sol2, 1.2)

Expand All @@ -120,7 +120,7 @@ get_normal_form(brpo_pd, 2)

# codim 2 PD
opts_pocoll_pd = ContinuationPar(brpo_pd.contparams, detect_bifurcation = 3, max_steps = 2, p_min = 1.e-2, dsmax = 1e-2, ds = 1e-3)
@set! opts_pocoll_pd.newton_options.tol = 1e-10
@reset opts_pocoll_pd.newton_options.tol = 1e-10
pd_po_coll2 = continuation(brpo_pd, 2, (@optic _.b0), opts_pocoll_pd;
verbosity = 0, plot = false,
detect_codim2_bifurcation = 1,
Expand Down
32 changes: 16 additions & 16 deletions test/codim2PO-shooting-mf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,15 @@ z0 = [0.1,0.1,1,0]
prob = BifurcationProblem(Pop, z0, par_pop, (@optic _.b0); record_from_solution = (x, p) -> (x = x[1], y = x[2], u = x[3]))

opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 4)
@set! opts_br.newton_options.verbose = true
@reset opts_br.newton_options.verbose = true

################################################################################
using OrdinaryDiffEq
prob_de = ODEProblem(Pop!, z0, (0, 600.), par_pop)
alg = Rodas5()
sol = solve(prob_de, alg)
sol = OrdinaryDiffEq.solve(prob_de, alg)
prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-10, abstol = 1e-12)
sol = solve(prob_de, Rodas5())
sol = OrdinaryDiffEq.solve(prob_de, Rodas5())
################################################################################
@info "plotting function"
argspo = (record_from_solution = (x, p; k...) -> begin
Expand Down Expand Up @@ -62,7 +62,7 @@ function flow(x0, prob0, tm, p = prob0.p)
end

@info "set AD"
# @set! probsh.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
# @reset probsh.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]

@info "Newton"
lspo = GMRESIterativeSolvers(verbose = false, N = length(cish), abstol = 1e-12, reltol = 1e-10)
Expand All @@ -76,7 +76,7 @@ _sol = BK.get_periodic_orbit(probsh, solpo.u, sol.prob.p)

@info "PO cont1"
opts_po_cont = setproperties(opts_br, max_steps = 50, save_eigenvectors = true, detect_loop = true, tol_stability = 1e-3, newton_options = optnpo)
@set! opts_po_cont.newton_options.verbose = false
@reset opts_po_cont.newton_options.verbose = false
br_fold_sh = continuation(probsh, cish, PALC(tangent = Bordered()), opts_po_cont;
# verbosity = 3, plot = true,
linear_algo = MatrixFreeBLS(lspo),
Expand All @@ -96,10 +96,10 @@ brpo_pd_sh = continuation(probsh2, cish, PALC(), opts_po_cont;
# codim 2 Fold
@info "--> Fold curve"
opts_posh_fold = ContinuationPar(br_fold_sh.contparams, detect_bifurcation = 0, max_steps = 0, p_min = 0.01, p_max = 1.2)
@set! opts_posh_fold.newton_options.tol = 1e-8
# @set! opts_posh_fold.newton_options.linsolver.solver.N = opts_posh_fold.newton_options.linsolver.solver.N+1
@set! opts_posh_fold.newton_options.verbose = false
@set! opts_posh_fold.newton_options.linsolver.solver.verbose=0
@reset opts_posh_fold.newton_options.tol = 1e-8
# @reset opts_posh_fold.newton_options.linsolver.solver.N = opts_posh_fold.newton_options.linsolver.solver.N+1
@reset opts_posh_fold.newton_options.verbose = false
@reset opts_posh_fold.newton_options.linsolver.solver.verbose=0
fold_po_sh1 = continuation(br_fold_sh, 2, (@optic _.ϵ), opts_posh_fold;
# verbosity = 3, #plot = true,
detect_codim2_bifurcation = 0,
Expand All @@ -114,7 +114,7 @@ fold_po_sh1 = continuation(br_fold_sh, 2, (@optic _.ϵ), opts_posh_fold;
# codim 2 PD
@info "--> PD curve"
opts_posh_pd = ContinuationPar(brpo_pd_sh.contparams, detect_bifurcation = 0, max_steps = 4, p_min = -1.)
@set! opts_posh_pd.newton_options.tol = 1e-8
@reset opts_posh_pd.newton_options.tol = 1e-8
pd_po_sh = continuation(brpo_pd_sh, 1, (@optic _.b0), opts_posh_pd;
verbosity = 0, #plot = true,
detect_codim2_bifurcation = 0,
Expand All @@ -136,15 +136,15 @@ pd_po_sh = continuation(brpo_pd_sh, 1, (@optic _.b0), opts_posh_pd;
#####
# find the PD NS case
par_pop2 = @set par_pop.b0 = 0.45
sol2 = solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5())
# plot(sol2, xlims= (8,10))

probshns, ci = generate_ci_problem( ShootingProblem(M=3), re_make(prob, params = sol2.prob.p), remake(prob_de, p = par_pop2), sol2, 1.; alg = Rodas5(),
jacobian = BK.AutoDiffMF()
)

# @set! probshns.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]
# @reset probshns.flow.vjp = (x,p,dx,tm) -> AD.pullback_function(AD.ZygoteBackend(), z->flow(z, prob_de,tm,p), x)(dx)[1]

brpo_ns = continuation(probshns, ci, PALC(), ContinuationPar(opts_po_cont; max_steps = 50, ds = -0.001);
# verbosity = 3, plot = true,
Expand All @@ -158,9 +158,9 @@ brpo_ns = continuation(probshns, ci, PALC(), ContinuationPar(opts_po_cont; max_s
# codim 2 NS
@info "--> NS curve"
opts_posh_ns = ContinuationPar(brpo_ns.contparams, detect_bifurcation = 0, max_steps = 0, p_min = -0., p_max = 1.2)
@set! opts_posh_ns.newton_options.tol = 1e-8
@set! opts_posh_ns.newton_options.linsolver.solver.verbose = 0
@set! opts_posh_ns.newton_options.verbose = false
@reset opts_posh_ns.newton_options.tol = 1e-8
@reset opts_posh_ns.newton_options.linsolver.solver.verbose = 0
@reset opts_posh_ns.newton_options.verbose = false
ns_po_sh = continuation(brpo_ns, 1, (@optic _.ϵ), opts_posh_ns;
# verbosity = 2, plot = true,
detect_codim2_bifurcation = 0,
Expand Down
12 changes: 6 additions & 6 deletions test/codim2PO-shooting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,9 @@ opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_
using OrdinaryDiffEq
prob_de = ODEProblem(Pop!, z0, (0,600.), par_pop)
alg = Rodas5()
sol = solve(prob_de, alg)
sol = OrdinaryDiffEq.solve(prob_de, alg)
prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-8, abstol = 1e-10)
sol = solve(prob_de, Rodas5())
sol = OrdinaryDiffEq.solve(prob_de, Rodas5())
################################################################################
argspo = (record_from_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, p.p)
Expand Down Expand Up @@ -112,8 +112,8 @@ BK.isinplace(_pdma)
#####
# find the PD NS case
par_pop2 = @set par_pop.b0 = 0.45
sol2 = solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5())
# plot(sol2, xlims= (8,10))

probshns, ci = generate_ci_problem(ShootingProblem(M=3), re_make(prob, params = sol2.prob.p), remake(prob_de, p = par_pop2), sol2, 1.; alg = Rodas5())
Expand Down Expand Up @@ -175,8 +175,8 @@ BK.NSMALinearSolver(_solpo, _p1[1], _p1[2], _probns.prob, _param, copy(_x.u), 1.

# find the PD case
par_pop2 = @set par_pop.b0 = 0.45
sol2 = solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5())
sol2 = solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0, 1000)), Rodas5())
sol2 = OrdinaryDiffEq.solve(remake(sol2.prob, tspan = (0, 10), u0 = sol2[end]), Rodas5())
# plot(sol2, xlims= (8,10))

probshpd, ci = generate_ci_problem(ShootingProblem(M=3), re_make(prob, params = sol2.prob.p), remake(prob_de, p = par_pop2), sol2, 1.; alg = Rodas5())
Expand Down
6 changes: 3 additions & 3 deletions test/event.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# using Revise
# using Plots
using Test
using BifurcationKit, Setfield
using BifurcationKit

const BK = BifurcationKit
###################################################################################################
Expand Down Expand Up @@ -216,7 +216,7 @@ ev3 = BK.BifDetectEvent
eves1 = BK.SetOfEvents(ev1, ev2, ev3)

args2 = @set args[end].detect_event = 2
@set! kwargs.verbosity = 0
@reset kwargs.verbosity = 0
br = continuation(args2...; kwargs...,
event = eves1,
)
Expand All @@ -241,7 +241,7 @@ br = continuation(args2...; kwargs...,
testBranch(br)

eves3 = SetOfEvents(eves1, eves2)
@set! kwargs.verbosity = 0
@reset kwargs.verbosity = 0
br = continuation(args2...; kwargs...,
event = eves3,
)
Expand Down
8 changes: 4 additions & 4 deletions test/poincareMap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ u0 = [.001, .001]
prob = ODEProblem(Fsl!, u0, (0., 100.), par_sl)
algsl = KenCarp4()#Rodas4P()
####################################################################################################
sol = solve(prob, algsl, abstol =1e-9, reltol=1e-6)
sol = OrdinaryDiffEq.solve(prob, algsl, abstol =1e-9, reltol=1e-6)

function flowTS(x, t, pb; alg = algsl, kwargs...)
_pb = remake(pb; u0 = x, tspan = (zero(eltype(t)), t) )
Expand All @@ -51,7 +51,7 @@ centers = [zeros(2)]
sectionH(x, c, n) = dot( x .- c[1], n[1])
pSection(u, t, integrator) = sectionH(u, centers, normals) * (integrator.iter > 1)
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(pSection, affect!; affect_neg! = nothing)
cb = OrdinaryDiffEq.ContinuousCallback(pSection, affect!; affect_neg! = nothing)

Π(x) = flowDE(x, Inf; callback = cb, save_everystep = false) # Poincaré return map
T(x) = flowTS(x, Inf64, prob; callback = cb)[1][end] # time to section
Expand Down Expand Up @@ -79,7 +79,7 @@ u0 = [0, 1.]
du0 = [0, -1.]

# check that we cross the sections the way we want
ts, ss = flowTS([0., 1], Inf, prob; callback = cb, save_everystep = true, saveat = LinRange(0,1,20))
ts, ss = flowTS([0., 1], Inf, prob; callback = cb, save_everystep = true, saveat = LinRange(0,1,20));
# plot(ts,ss[:,:]')
# plot(ss[1,:], ss[2,:], label="flow");scatter!(ss[1,[1]], ss[2,[1]]);plot!(sol[1,:], sol[2,:], label="sol")

Expand All @@ -104,7 +104,7 @@ println("\n──> Norm of the difference = ", resAna - resFD |> norminf)
const FD = ForwardDiff
const BK = BifurcationKit

tΣ, solΣ = flowTS(u0, Inf64, prob; callback = cb)
tΣ, solΣ = flowTS(u0, Inf64, prob; callback = cb);
= tΣ[end]; solΣ = solΣ[end]
= FD.jacobian( x -> flowTS(x, tΣ, prob)[2][end], (u0))
F = Fsl(Π(u0), par_sl)
Expand Down
Loading

0 comments on commit 5b957b2

Please sign in to comment.