Skip to content

Commit

Permalink
correct tests
Browse files Browse the repository at this point in the history
  • Loading branch information
rveltz committed Nov 28, 2024
1 parent 3cb63ee commit 39af849
Show file tree
Hide file tree
Showing 8 changed files with 49 additions and 27 deletions.
1 change: 1 addition & 0 deletions src/Continuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,7 @@ setparam(it::ContIterable{Tkind, Tprob, Talg, T}, p0::T) where {Tkind, Tprob, Ta
# getters
@inline getlens(it::ContIterable) = getlens(it.prob)
@inline getalg(it::ContIterable) = it.alg
@inline getprob(it::ContIterable) = it.prob
@inline callback(it::ContIterable) = it.callback_newton
record_from_solution(it::ContIterable) = record_from_solution(it.prob)
plot_solution(it::ContIterable) = plot_solution(it.prob)
Expand Down
4 changes: 4 additions & 0 deletions src/Newton.jl
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,10 @@ struct cbMaxNorm{T}
end
(cb::cbMaxNorm)(state; k...) = (return state.residual < cb.maxres)

"""
cb = cbMaxNormAndΔp(maxres, δp)
Create a callback used to reject residuals larger than `cb.maxres` or parameter step larger than `δp` in the Newton iterations. See docs for [`newton`](@ref).
"""
struct cbMaxNormAndΔp{T}
maxres::T
δp::T
Expand Down
2 changes: 1 addition & 1 deletion src/continuation/Contbase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ function _step_size_control!(state, contparams::ContinuationPar, verbosity)
ds = state.ds
if converged(state) == false
if abs(ds) <= contparams.dsmin
@error "Failure to converge with given tolerance = $(contparams.newton_options.tol).\n You can decrease the tolerance or pass a different norm `normC`.\n We reached the smallest value [dsmin] valid for ds, namely $(contparams.dsmin).\n Stopping continuation at continuation step $(state.step)."
@error "Failure to converge with given tolerance = $(contparams.newton_options.tol).\n You can decrease the tolerance or pass a different norm using the argument `normC`.\n We reached the smallest value [dsmin] valid for ds, namely $(contparams.dsmin).\n Stopping continuation at continuation step $(state.step)."
# we stop the continuation
state.stopcontinuation = true
return
Expand Down
44 changes: 28 additions & 16 deletions src/periodicorbit/PeriodicOrbitCollocation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -856,7 +856,9 @@ end
get_periodic_orbit(prob::PeriodicOrbitOCollProblem, x, p::Real) = get_periodic_orbit(prob, x, setparam(prob, p))

# same function as above but for coping with mesh adaptation
@views function get_periodic_orbit(prob::PeriodicOrbitOCollProblem, x::NamedTuple{(:mesh, :sol, :_mesh, :ϕ ), Tuple{Vector{Tp}, Vector{Tp}, Vector{Tp}, Vector{Tp}}}, p::Real) where Tp
@views function get_periodic_orbit(prob::PeriodicOrbitOCollProblem,
x::NamedTuple{(:mesh, :sol, :_mesh, :ϕ ), Tuple{Vector{Tp}, Vector{Tp}, Vector{Tp}, Vector{Tp}}},
p::Real) where Tp
mesh = x.mesh
u = x.sol
T = getperiod(prob, u, p)
Expand All @@ -865,7 +867,14 @@ get_periodic_orbit(prob::PeriodicOrbitOCollProblem, x, p::Real) = get_periodic_o
end

# function needed for automatic Branch switching from Hopf bifurcation point
function re_make(prob::PeriodicOrbitOCollProblem, prob_vf, hopfpt, ζr::AbstractVector, orbitguess_a, period; orbit = t -> t, k...)
function re_make(prob::PeriodicOrbitOCollProblem,
prob_vf,
hopfpt,
ζr::AbstractVector,
orbitguess_a,
period;
orbit = t -> t,
k...)
M = length(orbitguess_a)
N = length(ζr)

Expand Down Expand Up @@ -910,10 +919,10 @@ const DocStrjacobianPOColl = """
"""

function _newton_pocoll(probPO::PeriodicOrbitOCollProblem,
orbitguess,
options::NewtonPar;
defOp::Union{Nothing, DeflationOperator{T, Tf, vectype}} = nothing,
kwargs...) where {T, Tf, vectype}
orbitguess,
options::NewtonPar;
defOp::Union{Nothing, DeflationOperator{T, Tf, vectype}} = nothing,
kwargs...) where {T, Tf, vectype}
jacobianPO = probPO.jacobian
@assert jacobianPO in
(AutoDiffDense(), DenseAnalytical(), FullSparse(), DenseAnalyticalInplace()) "This jacobian $jacobianPO is not defined. Please chose another one."
Expand Down Expand Up @@ -980,7 +989,10 @@ newton(probPO::PeriodicOrbitOCollProblem,
_newton_pocoll(probPO, orbitguess, options; defOp = defOp, kwargs...)


function build_jacobian(coll::PeriodicOrbitOCollProblem, orbitguess, par; δ = convert(eltype(orbitguess), 1e-8))
function build_jacobian(coll::PeriodicOrbitOCollProblem,
orbitguess,
par;
δ = convert(eltype(orbitguess), 1e-8))
jacobianPO = coll.jacobian
@assert jacobianPO in (AutoDiffDense(), DenseAnalytical(), FullSparse(), FullSparseInplace(), DenseAnalyticalInplace()) "This jacobian is not defined. Please chose another one."

Expand Down Expand Up @@ -1062,10 +1074,10 @@ function continuation(probPO::PeriodicOrbitOCollProblem,
probwp = WrapPOColl(probPO, jacPO, orbitguess, getparams(probPO), getlens(probPO), _plotsol, _recordsol)

br = continuation(probwp, alg,
contParams;
kwargs...,
kind = PeriodicOrbitCont(),
finalise_solution = _finsol)
contParams;
kwargs...,
kind = PeriodicOrbitCont(),
finalise_solution = _finsol)
return br
end

Expand Down Expand Up @@ -1137,11 +1149,11 @@ References:
[2] R. D. Russell and J. Christiansen, “Adaptive Mesh Selection Strategies for Solving Boundary Value Problems,” SIAM Journal on Numerical Analysis 15, no. 1 (February 1978): 59–80, https://doi.org/10.1137/0715004.
"""
function compute_error!(coll::PeriodicOrbitOCollProblem, x::AbstractVector{Ty};
normE = norminf,
verbosity::Bool = false,
K = Inf,
par = nothing,
kw...) where Ty
normE = norminf,
verbosity::Bool = false,
K = Inf,
par = nothing,
kw...) where Ty
n, m, Ntst = size(coll) # recall that m = ncol
period = getperiod(coll, x, nothing)
# get solution, we copy x because it is overwritten at the end of this function
Expand Down
14 changes: 9 additions & 5 deletions test/codim2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -131,15 +131,19 @@ hp = newton(br, 2;
hp = newton(br, 2; options = NewtonPar( opts_br.newton_options; max_iterations = 10),start_with_eigen=true)

for eigen_start in (true, false)
hp_br = continuation(br, 2, (@optic _.k), ContinuationPar(opts_br, ds = -0.001, p_max = 1., p_min = 0., detect_bifurcation = 1, max_steps = 50, save_sol_every_step = 1, detect_event = 2), bdlinsolver = MatrixBLS(), start_with_eigen = eigen_start, update_minaug_every_step = 1, verbosity=0, plot=false)
hp_br = continuation(br, 2, (@optic _.k),
ContinuationPar(opts_br, ds = -0.001, p_max = 1., p_min = 0., detect_bifurcation = 1, max_steps = 50, save_sol_every_step = 1, detect_event = 2), bdlinsolver = MatrixBLS(),
start_with_eigen = eigen_start,
update_minaug_every_step = 1,
verbosity = 0,
detect_codim2_bifurcation = 2,
plot=false)
@test hp_br.kind isa BK.HopfCont
@test hp_br.specialpoint[1].type == :gh
@test hp_br.specialpoint[2].type == :nd
@test hp_br.specialpoint[3].type == :gh
@test hp_br.specialpoint[2].type == :gh

@test hp_br.specialpoint[1].param 0.305873681159479 rtol = 1e-5
@test hp_br.specialpoint[2].param 0.16452182436723148 atol = 1e-2
@test hp_br.specialpoint[3].param 0.23255761094689315 atol = 1e-4
@test hp_br.specialpoint[2].param 0.23255761094689315 atol = 1e-4

@test ~isnothing(hp_br.eig)
end
6 changes: 3 additions & 3 deletions test/codim2PO-shooting-mf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,9 @@ argspo = (record_from_solution = (x, p; k...) -> begin
period = getperiod(p.prob, x, p.p))
end,)
################################################################################
@info "import AD"
# @info "import AD"
# import AbstractDifferentiation as AD
@info "import Zygote"
# @info "import Zygote"
# using Zygote, SciMLSensitivity

@info "generate shooting problem"
Expand All @@ -61,7 +61,7 @@ function flow(x0, prob0, tm, p = prob0.p)
return sol[end]
end

@info "set AD"
# @info "set AD"
# @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"
Expand Down
3 changes: 2 additions & 1 deletion test/simple_continuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ x0 = 0.01 * ones(N)

# test continuation interface
begin
prob = BK.BifurcationProblem(F_simple, zeros(10), -1.5, (@optic _); J = Jac_simple)
iter = ContIterable(prob, PALC(), opts)
state = iterate(iter)
BK.gettangent(state[1])
Expand All @@ -101,7 +102,7 @@ end
###############
# basic continuation, without saving much information
_prob0 = BK.BifurcationProblem(F0_simple, [0.], -1.5, (@optic _))
opts0 = ContinuationPar(detect_bifurcation = 0, p_min = -2., save_sol_every_step = 0, max_steps = 40)
opts0 = ContinuationPar(detect_bifurcation = 0, p_min = -2., save_sol_every_step = 0, max_steps = 40, nev = 1)
_br0 = @time continuation(_prob0, PALC(), opts0) #597 allocations: 38.547 KiB

@test BK._hasstability(_br0) == false
Expand Down
2 changes: 1 addition & 1 deletion test/stuartLandauSH.jl
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ probPsh(initpo_bar, par_hopf)
# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finite_differences( x -> probPsh(x, par_hopf), initpo_bar)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-3
@test norm(_Jad - _Jana, Inf) < 3e-3

@info "Test newton"
ls = DefaultLS()
Expand Down

0 comments on commit 39af849

Please sign in to comment.