update readme

README.md

@@ -15,9 +15,9 @@ It incorporates continuation algorithms (PALC, deflated continuation, ...) based
 
 > The idea is to be able to seemingly switch the continuation algorithm a bit like changing the time stepper (Euler, RK4,...) for ODEs.
 
-`BifurcationKit` can also seek for periodic orbits of Cauchy problems. **It is by now, one of the only softwares which provides shooting methods *and* methods based on finite differences or collocation to compute periodic orbits.**
+`BifurcationKit` can also seek for [periodic orbits](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/) of Cauchy problems. **It is by now, one of the only software which provides shooting methods *and* methods based on finite differences / collocation to compute periodic orbits.**
 
-The current focus is on large scale nonlinear problems and multiple hardwares. Hence, the goal is to use Matrix Free methods on **GPU** (see [PDE example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2b/#The-Swift-Hohenberg-equation-on-the-GPU-(non-local)-1) and [Periodic orbit example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialsCGL/#Continuation-of-periodic-orbits-on-the-GPU-(Advanced)-1)) or on a **cluster** to solve non linear PDE, nonlocal problems, compute sub-manifolds...
+The current focus is on large scale nonlinear problems and multiple hardwares. Hence, the goal is to provide Matrix Free methods on **GPU** (see [PDE example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2b/#The-Swift-Hohenberg-equation-on-the-GPU-(non-local)-1) and [Periodic orbit example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialsCGL/#Continuation-of-periodic-orbits-on-the-GPU-(Advanced)-1)) or on **cluster** to study non linear PDE, nonlocal problems, compute sub-manifolds...
 
 > Despite this focus, the package can easily handle low dimensional problems and specific optimizations are regularly added.
 
@@ -56,7 +56,7 @@ To install the bleeding edge version, please run
 
 Most plugins are located in the organization [bifurcationkit](https://github.com/bifurcationkit):
 
-- [HclinicBifurcationKit.jl](https://github.com/bifurcationkit/HclinicBifurcationKit.jl) computation and bifurcation analysis of homoclinic / heteroclinic orbits of ordinary differential equations (ODE)
+- [HclinicBifurcationKit.jl](https://github.com/bifurcationkit/HclinicBifurcationKit.jl) bifurcation analysis of homoclinic / heteroclinic orbits of ordinary differential equations (ODE)
 - [DDEBifurcationKit.jl](https://github.com/bifurcationkit/DDEBifurcationKit.jl) bifurcation analysis of delay differential equations (DDE)
 - [AsymptoticNumericalMethod.jl](https://github.com/bifurcationkit/AsymptoticNumericalMethod.jl) provides the numerical continuation algorithm **Asymptotic Numerical Method** (ANM) which can be used directly in `BifurcationKit.jl`
 - [GridapBifurcationKit.jl](https://github.com/bifurcationkit/GridapBifurcationKit) bifurcation analysis of PDEs solved with the Finite Elements Method (FEM) using the package [Gridap.jl](https://github.com/gridap/Gridap.jl).

src/continuation/MoorePenrose.jl

@@ -32,8 +32,8 @@ internal_adaptation!(alg::MoorePenrose, swch::Bool) = internal_adaptation!(alg.t
 """
 $(SIGNATURES)
 """
-function MoorePenrose(;tangent = PALC(), 
-                        method = direct, 
+function MoorePenrose(;tangent = PALC(),
+                        method = direct,
                         ls = nothing)
     if ~(method == iterative)
         ls = isnothing(ls) ? DefaultLS() : ls
@@ -120,7 +120,7 @@ function newton_moore_penrose(iter::AbstractContinuationIterable,
     ϵ = getdelta(prob)
     paramlens = getlens(iter)
     contparams = getcontparams(iter)
-    T = eltype(iter)
+    𝒯 = eltype(iter)
 
     @unpack method = iter.alg
 
@@ -142,7 +142,7 @@ function newton_moore_penrose(iter::AbstractContinuationIterable,
     dX = _copy(res_f) # copy(res_f)
     # dFdp = (F(x, p + ϵ) - res_f) / ϵ
     dFdp = _copy(residual(prob, x, set(par, paramlens, p + ϵ)))
-    minus!(dFdp, res_f); rmul!(dFdp, T(1) / ϵ)
+    minus!(dFdp, res_f); rmul!(dFdp, one(𝒯) / ϵ)
 
     res = normN(res_f)
     residuals = [res]
@@ -162,35 +162,37 @@ function newton_moore_penrose(iter::AbstractContinuationIterable,
     X = BorderedArray(x, p)
     if linsolver isa AbstractIterativeLinearSolver || (method == iterative)
         ϕ = _copy(τ0)
-        rmul!(ϕ,  T(1) / norm(ϕ))
+        rmul!(ϕ,  one(𝒯) / norm(ϕ))
     end
 
     while (step < max_iterations) && (res > tol) && line_step && compute
         step += 1
         # dFdp = (F(x, p + ϵ) - F(x, p)) / ϵ)
         copyto!(dFdp, residual(prob, x, set(par, paramlens, p + ϵ)))
-        minus!(dFdp, res_f); rmul!(dFdp, T(1) / ϵ)
+        minus!(dFdp, res_f); rmul!(dFdp, one(𝒯) / ϵ)
 
         # compute jacobian
         J = jacobian(prob, x, set(par, paramlens, p))
         if method == direct || method == pInv
             @debug "Moore-Penrose direct/pInv"
             Jb = hcat(J, dFdp)
             if method == direct
-                dx, flag, converged = linsolver(Jb, res_f)
+                dx, flag, itlinear = linsolver(Jb, res_f)
+                ~flag && @debug "[MoorePenrose] Linear solver did not converge."
             else
-                # pinv(Array(Jb)) * res_f seems to work better than the following
-                dx = LinearAlgebra.pinv(Array(Jb)) * res_f; flag = true;
+                # dx = pinv(Array(Jb)) * res_f #seems to work better than the following
+                dx = LinearAlgebra.pinv(Array(Jb)) * res_f
+                flag = true;
+                itlinear = 1
             end
             x .-= @view dx[begin:end-1]
             p -= dx[end]
-            itlinear = 1
         else
             @debug "Moore-Penrose Iterative"
             # A = hcat(J, dFdp); A = vcat(A, ϕ')
             # X .= X .- A \ vcat(res_f, 0)
             # x .= X[begin:end-1]; p = X[end]
-            du, dup, flag, itlinear1 = linsolver(J, dFdp, ϕ.u, ϕ.p, res_f, zero(T), one(T), one(T))
+            du, dup, flag, itlinear1 = linsolver(J, dFdp, ϕ.u, ϕ.p, res_f, zero(𝒯), one(𝒯), one(𝒯))
             minus!(x, du)
             p -= dup
             verbose && print_nonlinear_step(step, nothing, itlinear1)
@@ -204,13 +206,13 @@ function newton_moore_penrose(iter::AbstractContinuationIterable,
             # compute jacobian
             J = jacobian(prob, x, set(par, paramlens, p))
             copyto!(dFdp, residual(prob, x, set(par, paramlens, p + ϵ)))
-            minus!(dFdp, res_f); rmul!(dFdp, T(1) / ϵ)
+            minus!(dFdp, res_f); rmul!(dFdp, one(𝒯) / ϵ)
             # A = hcat(J, dFdp); A = vcat(A, ϕ')
             # ϕ .= A \ vcat(zero(x),1)
-            u, up, flag, itlinear2 = linsolver(J, dFdp, ϕ.u, ϕ.p, zero(x), one(T), one(T), one(T))
-            ~flag && @debug "Linear solver for (J-iω) did not converge."
+            u, up, flag, itlinear2 = linsolver(J, dFdp, ϕ.u, ϕ.p, zero(x), one(𝒯), one(𝒯), one(𝒯))
+            ~flag && @debug "[MoorePenrose] Linear solver did not converge."
             ϕ.u .= u; ϕ.p = up
-            # rmul!(ϕ,  T(1) / norm(ϕ))
+            # rmul!(ϕ,  one(𝒯) / norm(ϕ))
             itlinear = (itlinear1 .+ itlinear2)
         end
         push!(residuals, res)
@@ -222,5 +224,11 @@ function newton_moore_penrose(iter::AbstractContinuationIterable,
     end
     verbose && print_nonlinear_step(step, res, 0, true) # display last line of the table
     flag = (residuals[end] < tol) & callback((;x, res_f, nothing, residual=res, step, contparams, p, residuals, z0); fromNewton = false, kwargs...)
-    return NonLinearSolution(BorderedArray(x, p), prob, residuals, flag, step, itlineartot)
+
+    return NonLinearSolution(BorderedArray(x, p),
+                            prob,
+                            residuals,
+                            flag,
+                            step,
+                            itlineartot)
 end

src/continuation/Palc.jl

@@ -506,6 +506,12 @@ function newton_palc(iter::AbstractContinuationIterable,
         compute = callback((;x, res_f, J, residual=res, step, itlinear, contparams, z0, p, residuals, options = (;linsolver)); fromNewton = false, kwargs...)
     end
     verbose && print_nonlinear_step(step, res, 0, true) # display last line of the table
-    flag = (residuals[end] < tol) & callback((;x, res_f, residual=res, step, contparams, p, residuals, options = (;linsolver)); fromNewton = false, kwargs...)
-    return NonLinearSolution(BorderedArray(x, p), prob, residuals, flag, step, itlineartot)
+    flag = (residuals[end] < tol) & callback((;x, res_f, residual = res, step, contparams, p, residuals, options = (;linsolver)); fromNewton = false, kwargs...)
+
+    return NonLinearSolution(BorderedArray(x, p),
+                            prob,
+                            residuals,
+                            flag,
+                            step,
+                            itlineartot)
 end

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -818,7 +818,7 @@ 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)
 
@@ -1006,21 +1006,23 @@ function getmaximum(prob::PeriodicOrbitOCollProblem, x::AbstractVector, p)
 end
 
 # this function updates the section during the continuation run
-@views function updatesection!(prob::PeriodicOrbitOCollProblem, x, par; stride = 0)
+@views function updatesection!(prob::PeriodicOrbitOCollProblem, 
+                                x::AbstractVector, 
+                                par)
     @debug "Update section Collocation"
     # update the reference point
     prob.xπ .= 0
 
     # update the "normals"
-    prob.ϕ .= x[1:end-1]
+    prob.ϕ .= x[begin:end-1]
     return true
 end
 ####################################################################################################
 # mesh adaptation method
 
 # iterated derivatives
 ∂(f) = x -> ForwardDiff.derivative(f, x)
-∂(f, n) = n == 0 ? f : ∂(∂(f), n-1)
+∂(f, n::Int) = n == 0 ? f : ∂(∂(f), n-1)
 
 @views function (sol::POSolution{ <: PeriodicOrbitOCollProblem})(t0)
     n, m, Ntst = size(sol.pb)