Skip to content

Commit

Permalink
use the field in_bisection of ContState. Remove it from kwargs to fin…
Browse files Browse the repository at this point in the history
…alizer
  • Loading branch information
rveltz committed Feb 4, 2024
1 parent 196fef8 commit b94e569
Show file tree
Hide file tree
Showing 12 changed files with 36 additions and 28 deletions.
2 changes: 1 addition & 1 deletion src/Bifurcations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@ function locate_bifurcation!(iter::ContIterable, _state::ContState, verbose::Boo
state.step > 0 && (interval = @set interval[indinterval] = getp(state))

# we call the finalizer
# iter.finalise_solution(state.z_old, state.tau, state.step, nothing; bisection = true)
# iter.finalise_solution(state.z_old, state.tau, state.step, nothing; )

if verbose
ct0 = rightmost(state.eigvals)
Expand Down
2 changes: 2 additions & 0 deletions src/Continuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,7 @@ function Base.copyto!(dest::ContState, src::ContState)
end

# getters
@inline converged(::Nothing) = false
@inline converged(state::AbstractContinuationState) = state.converged
getsolution(state::AbstractContinuationState) = state.z
getx(state::AbstractContinuationState) = state.z.u
Expand All @@ -202,6 +203,7 @@ getx(state::AbstractContinuationState) = state.z.u
@inline is_stable(state::AbstractContinuationState) = state.n_unstable[1] == 0
@inline stepsizecontrol(state::AbstractContinuationState) = state.stepsizecontrol
@inline in_bisection(state::AbstractContinuationState) = state.in_bisection
@inline in_bisection(::Nothing) = false
####################################################################################################
# condition for halting the continuation procedure (i.e. when returning false)
@inline done(it::ContIterable, state::ContState) =
Expand Down
4 changes: 2 additions & 2 deletions src/DeflatedContinuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -410,13 +410,13 @@ end
# _ps = getParamValues(brs)
# display(_ps)
# if δp > 0
# @assert 1==1 "searchsorted?"
# @assert true "searchsorted?"
# idp = findfirst(x -> x.p > startingpoint.p, _ps)
# @show idp
# pspan = (startingpoint.p, ~isnothing(idp) ? _ps[idp].p : pMax)
# printstyled(color=:blue, "──▶ to the right\n")
# else
# @assert 1==1 "searchsorted?"
# @assert true "searchsorted?"
# idp = findfirst(x -> x.p < startingpoint.p, _ps)
# pspan = (~isnothing(idp) ? _ps[idp].p : pMin, startingpoint.p)
# printstyled(color=:green, "──▶ to the left\n")
Expand Down
2 changes: 1 addition & 1 deletion src/events/EventDetection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ function locate_event!(event::AbstractEvent, iter, _state, verbose::Bool = true)
state.step > 0 && (@set! interval[indinterval] = getp(state))

# we call the finalizer
state.stopcontinuation = ~iter.finalise_solution(state.z, state.τ, state.step, nothing; bisection = true, state = state)
state.stopcontinuation = ~iter.finalise_solution(state.z, state.τ, state.step, nothing; state = state)

if verbose
printstyled(color=:blue, bold = true, "────> ", state.step,
Expand Down
21 changes: 12 additions & 9 deletions src/periodicorbit/PeriodicOrbitUtils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ The arguments are
- `M` number of time slices in the periodic orbit guess
- `amplitude`: amplitude of the periodic orbit guess
"""
function guess_from_hopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)
function guess_from_hopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M::Int, amplitude; phase = 0)
hopfpoint = HopfPoint(br, ind_hopf)
specialpoint = br.specialpoint[ind_hopf]

Expand All @@ -29,7 +29,7 @@ function guess_from_hopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplit

# vec_hopf is the eigenvector for the eigenvalues iω
vec_hopf = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev-1)
vec_hopf ./= norm(vec_hopf)
vec_hopf ./= norm(vec_hopf)

orbitguess = [real.(hopfpoint.u .+ amplitude .* vec_hopf .* exp(-2pi * complex(0, 1) .* (ii/(M-1) - phase))) for ii in 0:M-1]

Expand All @@ -53,12 +53,14 @@ function modify_po_finalise(prob, kwargs, updateSectionEveryStep)
return Finaliser(prob, get(kwargs, :finalise_solution, nothing), updateSectionEveryStep)
end

function (finalizer::Finaliser{ <: AbstractPeriodicOrbitProblem})(z, tau, step, contResult; bisection = false, kF...)
function (finalizer::Finaliser{ <: AbstractPeriodicOrbitProblem})(z, tau, step, contResult; kF...)
updateSectionEveryStep = finalizer.updateSectionEveryStep
# we first check that the continuation step was successful
# if not, we do not update the problem with bad information
success = converged(get(kF, :state, nothing))
if success && mod_counter(step, updateSectionEveryStep) == 1 && ~bisection
state = get(kF, :state, nothing)
success = converged(state)
bisection = in_bisection(state)
if success && mod_counter(step, updateSectionEveryStep) == 1 && bisection == false
@debug "[Periodic orbit] update section"
# Trapezoid and Shooting need the parameters for section update:
updatesection!(finalizer.prob, z.u, setparam(contResult, z.p))
Expand All @@ -72,15 +74,16 @@ end

# version specific to collocation. Handle mesh adaptation
function (finalizer::Finaliser{ <: Union{ <: PeriodicOrbitOCollProblem,
<: WrapPOColl}})(z, tau, step, contResult; bisection = false, kF...)
<: WrapPOColl}})(z, tau, step, contResult; kF...)
updateSectionEveryStep = finalizer.updateSectionEveryStep
coll = finalizer.prob
# we first check that the continuation step was successful
# if not, we do not update the problem with bad information
state = get(kF, :state, nothing)
success = isnothing(state) ? false : converged(state)
success = converged(state)
bisection = in_bisection(state)
# mesh adaptation
if success && coll.meshadapt && ~bisection
if success && coll.meshadapt && bisection == false
@debug "[Collocation] update mesh"
oldsol = _copy(z) # avoid possible overwrite in compute_error!
oldmesh = get_times(coll) .* getperiod(coll, oldsol.u, nothing)
Expand All @@ -92,7 +95,7 @@ function (finalizer::Finaliser{ <: Union{ <: PeriodicOrbitOCollProblem,
return false
end
end
if success && mod_counter(step, updateSectionEveryStep) == 1 && ~bisection
if success && mod_counter(step, updateSectionEveryStep) == 1 && bisection == false
@debug "[collocation] update section"
updatesection!(coll, z.u, setparam(contResult, z.p))
end
Expand Down
2 changes: 1 addition & 1 deletion src/periodicorbit/PoincareRM.jl
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ function _solve(Π::PoincaréMap{ <: WrapPOSh}, xₛ, par)
x₀ = Π.po[N+1:end]
x₀[end] = sh.ds[end]
mapΠ(x, p) = poincaré_functional(Π, x, p, xₛ)
# @assert 1==0 "needs a jacobian"
# @assert false "needs a jacobian"
probΠ = BifurcationProblem(mapΠ,
x₀,
par)
Expand Down
2 changes: 1 addition & 1 deletion src/periodicorbit/PoincareShooting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -338,7 +338,7 @@ function (psh::PoincareShootingProblem)(x_bar::AbstractVector, par, dx_bar::Abst
@views outc[:, ii] .= dxc[:, ii] .- diff_poincare_map(psh, xc[:, im1], par, dxc[:, im1], im1)
end
else
@assert 1==0 "Analytical Jacobian for parallel Poincare Shooting not implemented yet. Please use the option δ > 0 to use Matrix-Free jacobian or chose `:FiniteDifferencesDense` to compute jacobian based on finite differences."
@assert false "Analytical Jacobian for parallel Poincare Shooting not implemented yet. Please use the option δ > 0 to use Matrix-Free jacobian or chose `:FiniteDifferencesDense` to compute jacobian based on finite differences."
end

# build the array to be returned
Expand Down
4 changes: 2 additions & 2 deletions src/periodicorbit/StandardShooting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ function (sh::ShootingProblem)(x::BorderedArray, par)
copyto!(out.u[ii], evolve(sh.flow, xc[ii], par, sh.ds[ii] * T).u .- xc[ip1])
end
else
@assert 1==0 "Not implemented yet. Try to use an AbstractVector instead"
@assert false "Not implemented yet. Try to use an AbstractVector instead"
end

# add constraint
Expand Down Expand Up @@ -245,7 +245,7 @@ function (sh::ShootingProblem)(x::BorderedArray, par, dx::BorderedArray; δ = co
copyto!(out.u[ii], tmp.du .+ vf(sh.flow, tmp.u, par) .* sh.ds[ii] .* dT .- dx.u[ip1])
end
else
@assert 1==0 "Not implemented yet. Try using AbstractVectors instead"
@assert false "Not implemented yet. Try using AbstractVectors instead"
end

# add constraint
Expand Down
4 changes: 2 additions & 2 deletions src/periodicorbit/codim2/MinAugNS.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ end

function apply_jacobian_neimark_sacker(pb, x, par, ω, dx, _transpose = false)
if _transpose == false
@assert 1==0
@assert false
return jacobian_adjoint_neimark_sacker_matrix_free(pb, x, par, ω, dx)
else
# if matrix-free:
Expand Down Expand Up @@ -197,7 +197,7 @@ function NSMALinearSolver(x, p::𝒯, ω::𝒯, 𝐍𝐒::NeimarkSackerProblemMi

return x1 .- dp .* x2, dp, dω, true, it1 + it2 + sum(itv) + sum(itw)
else
@assert 1==0 "WIP. Please select another jacobian method like :autodiff or :finiteDifferences. You can also pass the option usehessian = false."
@assert false "WIP. Please select another jacobian method like :autodiff or :finiteDifferences. You can also pass the option usehessian = false."
end

return dX, dsig, true, sum(it) + sum(itv) + sum(itw)
Expand Down
2 changes: 1 addition & 1 deletion src/periodicorbit/codim2/MinAugPD.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ function apply_jacobian_period_doubling(pb, x, par, dx, _transpose = false)
# else
# return apply(jacobian_period_doubling(pb, x, par), dx)
# end
@assert 1==0 "Please report to the website of BifurcationKit"
@assert false "Please report to the website of BifurcationKit"
else
# if matrix-free:
if has_adjoint(pb)
Expand Down
4 changes: 2 additions & 2 deletions src/periodicorbit/codim2/StandardShooting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ function jacobian_pd_nf_matrix_free(pbwrap::WrapPOSh{ <: ShootingProblem }, x, p
end
end
else
@assert 1==0 "WIP! No parallel matrix-free shooting for curve of PD/NS"
@assert false "WIP! No parallel matrix-free shooting for curve of PD/NS"
# call jacobian of the flow, jacobian-vector product
solOde = jvp(sh.flow, xc, par, dxc, sh.ds .* T)
for ii in 1:M
Expand Down Expand Up @@ -110,7 +110,7 @@ function jacobian_adjoint_pd_nf_matrix_free(pbwrap::WrapPOSh{ <: ShootingProblem
end
end
else
@assert 1==0
@assert false
# call jacobian of the flow, jacobian-vector product
solOde = jvp(sh.flow, xc, par, dxc, sh.ds .* T)
for ii in 1:M
Expand Down
15 changes: 9 additions & 6 deletions src/periodicorbit/codim2/codim2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -136,8 +136,10 @@ function (finalizer::Finaliser{<: AbstractMABifurcationProblem})(z, tau, step, c
updateSectionEveryStep = finalizer.updateSectionEveryStep
# we first check that the continuation step was successful
# if not, we do not update the problem with bad information
success = converged(get(kF, :state, nothing))
if success && mod_counter(step, updateSectionEveryStep) == 1 && ~bisection
state = get(kF, :state, nothing)
success = converged(state)
bisection = in_bisection(state)
if success && mod_counter(step, updateSectionEveryStep) == 1 && bisection == false
# we get the MA problem
wrap_ma = finalizer.prob
𝐏𝐛 = wrap_ma.prob
Expand All @@ -159,7 +161,7 @@ function (finalizer::Finaliser{<: AbstractMABifurcationProblem})(z, tau, step, c
end
end

function (finalizer::Finaliser{<: AbstractMABifurcationProblem{ <: AbstractProblemMinimallyAugmented{ <: WrapPOColl}}})(Z, tau, step, contResult; bisection = false, kF...)
function (finalizer::Finaliser{<: AbstractMABifurcationProblem{ <: AbstractProblemMinimallyAugmented{ <: WrapPOColl}}})(Z, tau, step, contResult; kF...)
updateSectionEveryStep = finalizer.updateSectionEveryStep
𝐏𝐛 = finalizer.prob.prob
coll = 𝐏𝐛.prob_vf.prob
Expand All @@ -168,9 +170,10 @@ function (finalizer::Finaliser{<: AbstractMABifurcationProblem{ <: AbstractProbl
# we first check that the continuation step was successful
# if not, we do not update the problem with bad information
state = get(kF, :state, nothing)
success = isnothing(state) ? false : converged(state)
success = converged(state)
bisection = in_bisection(state)
# mesh adaptation
if success && coll.meshadapt && ~bisection
if success && coll.meshadapt && bisection == false
@debug "[Collocation] update mesh"
oldsol = _copy(x) # avoid possible overwrite in compute_error!
oldmesh = get_times(coll) .* getperiod(coll, oldsol, nothing)
Expand All @@ -182,7 +185,7 @@ function (finalizer::Finaliser{<: AbstractMABifurcationProblem{ <: AbstractProbl
return false
end
end
if success && mod_counter(step, updateSectionEveryStep) == 1 && ~bisection
if success && mod_counter(step, updateSectionEveryStep) == 1 && bisection == false
@debug "[collocation] update section"
updatesection!(coll, x, nothing) # collocation does not need the parameter for updatesection!
end
Expand Down

0 comments on commit b94e569

Please sign in to comment.