From 39af849b96310168e9d68a3f4ea6c9e91c79bc69 Mon Sep 17 00:00:00 2001 From: romain veltz Date: Thu, 28 Nov 2024 22:41:15 +0100 Subject: [PATCH] correct tests --- src/Continuation.jl | 1 + src/Newton.jl | 4 ++ src/continuation/Contbase.jl | 2 +- src/periodicorbit/PeriodicOrbitCollocation.jl | 44 ++++++++++++------- test/codim2.jl | 14 +++--- test/codim2PO-shooting-mf.jl | 6 +-- test/simple_continuation.jl | 3 +- test/stuartLandauSH.jl | 2 +- 8 files changed, 49 insertions(+), 27 deletions(-) diff --git a/src/Continuation.jl b/src/Continuation.jl index 64e62c4c..a1f0b5b8 100644 --- a/src/Continuation.jl +++ b/src/Continuation.jl @@ -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) diff --git a/src/Newton.jl b/src/Newton.jl index 807461f8..6bc728b8 100644 --- a/src/Newton.jl +++ b/src/Newton.jl @@ -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 diff --git a/src/continuation/Contbase.jl b/src/continuation/Contbase.jl index e7c0e0a1..5b07e2b1 100644 --- a/src/continuation/Contbase.jl +++ b/src/continuation/Contbase.jl @@ -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 diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index cc83218b..3aa0c219 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/src/periodicorbit/PeriodicOrbitCollocation.jl @@ -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) @@ -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) @@ -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." @@ -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." @@ -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 @@ -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 diff --git a/test/codim2.jl b/test/codim2.jl index 1e58b6e1..5b208897 100644 --- a/test/codim2.jl +++ b/test/codim2.jl @@ -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 diff --git a/test/codim2PO-shooting-mf.jl b/test/codim2PO-shooting-mf.jl index fb75c94d..759137a9 100644 --- a/test/codim2PO-shooting-mf.jl +++ b/test/codim2PO-shooting-mf.jl @@ -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" @@ -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" diff --git a/test/simple_continuation.jl b/test/simple_continuation.jl index 1d0649e3..e3f0397f 100644 --- a/test/simple_continuation.jl +++ b/test/simple_continuation.jl @@ -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]) @@ -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 diff --git a/test/stuartLandauSH.jl b/test/stuartLandauSH.jl index 3975e12d..c1cc3b01 100644 --- a/test/stuartLandauSH.jl +++ b/test/stuartLandauSH.jl @@ -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()