From b3ba34c5815cc485c0fe12b2306a9cb9923786d3 Mon Sep 17 00:00:00 2001 From: romain veltz Date: Sat, 27 Jan 2024 15:01:21 +0100 Subject: [PATCH] update readme --- README.md | 6 +-- src/continuation/MoorePenrose.jl | 40 +++++++++++-------- src/continuation/Palc.jl | 10 ++++- src/periodicorbit/PeriodicOrbitCollocation.jl | 10 +++-- 4 files changed, 41 insertions(+), 25 deletions(-) diff --git a/README.md b/README.md index 10446d36..a4a2580c 100644 --- a/README.md +++ b/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). diff --git a/src/continuation/MoorePenrose.jl b/src/continuation/MoorePenrose.jl index af4fd45a..7589a77c 100644 --- a/src/continuation/MoorePenrose.jl +++ b/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,14 +162,14 @@ 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)) @@ -177,20 +177,22 @@ function newton_moore_penrose(iter::AbstractContinuationIterable, @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 diff --git a/src/continuation/Palc.jl b/src/continuation/Palc.jl index e3c8a6a5..748da48f 100644 --- a/src/continuation/Palc.jl +++ b/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 diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index e7458dd9..4b1ef5a0 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/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,13 +1006,15 @@ 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 #################################################################################################### @@ -1020,7 +1022,7 @@ end # 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)