diff --git a/README.md b/README.md index e15ad303..ec7918ba 100644 --- a/README.md +++ b/README.md @@ -62,11 +62,16 @@ Most plugins are located in the organization [bifurcationkit](https://github.com - [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). - [PeriodicSchurBifurcationKit.jl](https://github.com/bifurcationkit/PeriodicSchurBifurcationKit.jl) state of the art computation of Floquet coefficients, useful for computing the stability of periodic orbits. -## Examples of bifurcation diagrams +## Overview of capabilities +The list of capabilities is available [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/capabilities/). -| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/BDSH1d.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittlemannBD-1.png) | +## Examples of bifurcation diagrams (ODEs and PDEs) + +| ![](https://github.com/bifurcationkit/BifurcationKitDocs.jl/blob/main/docs/src/tutorials/ode/nm-per.png?raw=true) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/com-fig3.png?raw=true) | |:-------------:|:-------------:| +| [simple ODE example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialsODE/#Neural-mass-equation-(Hopf-aBS)) | [Codimension 2 (ODE)](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialCO/#CO-oxydation-(codim-2)) | +| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/BDSH1d.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittlemannBD-1.png) | | [Automatic Bif. Diagram in 1D Swift Hohenberg](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/Swift-Hohenberg1d/#d-Swift-Hohenberg-equation-(Automatic)) | [Automatic Bif. Diagram in 2D Bratu](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittelmannAuto/#Automatic-diagram-of-2d-Bratu–Gelfand-problem-(Intermediate)) | | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/sh2dbranches.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/bru-po-cont-3br.png) | | [Snaking in 2D Swift Hohenberg](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2/#d-Swift-Hohenberg-equation:-snaking,-Finite-Differences) | [Periodic orbits in 1D Brusselator](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials3/#d-Brusselator-(automatic)) @@ -75,7 +80,3 @@ Most plugins are located in the organization [bifurcationkit](https://github.com | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/carrier.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/GPU-branch.png) | | [Deflated Continuation in Carrier problem](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialCarrier/#Deflated-Continuation-in-the-Carrier-Problem) | [2D Swift Hohenberg on GPU](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2b/#d-Swift-Hohenberg-equation-(non-local)-on-the-GPU,-periodic-BC-(Advanced)) | - -## Main features - -The list of capabilities is available [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/capabilities/). diff --git a/examples/SHpde_snaking.jl b/examples/SHpde_snaking.jl index 1f13f1e5..1045aa86 100644 --- a/examples/SHpde_snaking.jl +++ b/examples/SHpde_snaking.jl @@ -37,9 +37,9 @@ prob = BifurcationProblem(R_SH, sol0, parSH, (@lens _.λ); J = Jac_sp, plotSolution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...))) #################################################################################################### optnew = NewtonPar(verbose = false, tol = 1e-12) - # allocations 357, 0.8ms - sol1 = @time newton(prob, optnew) - Plots.plot(X, sol1.u) +# allocations 357, 0.8ms +sol1 = @time newton(prob, optnew) +Plots.plot(X, sol1.u) opts = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds = 0.01, newtonOptions = setproperties(optnew; maxIter = 30, tol = 1e-8), pMax = 1., @@ -63,7 +63,7 @@ kwargsC = (verbosity = 3, callbackN = cb ) -brflat = @time continuation(prob, PALC(), opts; kwargsC...) +brflat = @time continuation(prob, PALC(), opts; kwargsC..., verbosity = 0) plot(brflat, putspecialptlegend = false) #################################################################################################### diff --git a/examples/mittleman.jl b/examples/mittleman.jl index ade87223..33bf619c 100644 --- a/examples/mittleman.jl +++ b/examples/mittleman.jl @@ -158,7 +158,7 @@ code = () title!("#branches = $(size(getBranch(diagram, code)))") # xlims!(0.01, 0.065, ylims=(2.5,6.5)) -plot(getBranchesFromBP(diagram, 2); plotfold = false, legend = false, vars = (:param, :x)) +plot(getBranchesFromBP(diagram, 2); plotfold = false, legend = false, vars = (:param, :n2)) getBranch(diagram, (2,1)) |> plot @@ -212,7 +212,6 @@ res = BK.continuation(br, 2, plot(res..., br ;plotfold= false) - δp = 0.005 deflationOp = DeflationOperator(2, 1.0, [zeros(2)]) success = [0] diff --git a/src/BifurcationKit.jl b/src/BifurcationKit.jl index 19e7cc3a..3c134d45 100644 --- a/src/BifurcationKit.jl +++ b/src/BifurcationKit.jl @@ -203,7 +203,7 @@ module BifurcationKit export NewtonPar, newton, newtonDeflated, newtonPALC, newtonFold, newtonHopf, newtonBordered, NonLinearSolution # continuation methods - export ContinuationPar, ContResult, GenericBifPoint, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered, eigenvec, eigenvals, getSolx, getSolp, bifurcation_points + export ContinuationPar, ContResult, GenericBifPoint, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered, eigenvec, eigenvals, getSolx, getSolp, bifurcation_points, SpecialPoint # events export ContinuousEvent, DiscreteEvent, PairOfEvents, SetOfEvents, SaveAtEvent, FoldDetectEvent, BifDetectEvent diff --git a/src/BifurcationPoints.jl b/src/BifurcationPoints.jl index b900da4c..de184b06 100644 --- a/src/BifurcationPoints.jl +++ b/src/BifurcationPoints.jl @@ -8,24 +8,34 @@ istranscritical(bp::AbstractBranchPoint) = false """ $(TYPEDEF) -Structure to record specials point on the curve. There are two types of specials point that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/). +Structure to record special points on a curve. There are two types of special points that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/). $(TYPEDFIELDS) """ @with_kw struct SpecialPoint{T, Tp, Tv, Tvτ} <: AbstractBifurcationPoint - "Description of the special point. In case of Events, this field records the user passed named to the event, or the default `:userD`, `:userC`. In case of bifurcation points, it can be the following: + "Description of the special points. In case of Events, this field records the user passed named to the event, or the default `:userD`, `:userC`. In case of bifurcation points, it can be one of the following: - :bp Bifurcation point, simple eigenvalue crossing the imaginary axis - :fold Fold point - :hopf Hopf point - :nd not documented bifurcation point. Detected by multiple eigenvalues crossing. Generally occurs in problems with symmetries or in cases where the continuation step size is too large and merge two different bifurcation points. - :cusp Cusp point - - :gh Generalized Hopf point + - :gh Generalized Hopf point (also called Bautin point) - :bt Bogdanov-Takens point - :zh Zero-Hopf point - :hh Hopf-Hopf point - :ns Neimark-Sacker point - :pd Period-doubling point + - :R1 Strong resonance 1:1 of periodic orbits + - :R2 Strong resonance 1:2 of periodic orbits + - :R3 Strong resonance 1:3 of periodic orbits + - :R4 Strong resonance 1:4 of periodic orbits + - :foldFlip Fold / Flip of periodic orbits + - :foldNS Fold / Neimark-Sacker of periodic orbits + - :pdNS Period-Doubling / Neimark-Sacker of periodic orbits + - :gpd Generalized Period-Doubling of periodic orbits + - :nsns Double Neimark-Sacker of periodic orbits + - :ch Chenciner bifurcation of periodic orbits " type::Symbol = :none @@ -62,7 +72,7 @@ $(TYPEDFIELDS) "Precision in the location of the special point" precision::T = -1 - "Interval containing the special point" + "Interval parameter containing the special point" interval::Tuple{T, T} = (0, 0) end @@ -98,7 +108,7 @@ SpecialPoint(it::ContIterable, state::ContState, type::Symbol, status::Symbol, i function _show(io::IO, bp::SpecialPoint, ii::Int, p::String = "p") if bp.type == :none ; return; end if bp.type == :endpoint - @printf(io, "- #%3i, \033[1m%5s\033[0m at %s ≈ %+4.8f, step = %3i\n", ii, "endpoint", p, bp.param, bp.step) + @printf(io, "- #%3i, \033[1m%5s\033[0m at %s ≈ %+4.8f, step = %3i\n", ii, "endpoint", p, bp.param, bp.step) return end if bp.status == :converged diff --git a/src/Problems.jl b/src/Problems.jl index 7c1eb909..ca30a6e3 100644 --- a/src/Problems.jl +++ b/src/Problems.jl @@ -17,14 +17,13 @@ $(TYPEDEF) Structure to hold the vector field and its derivatives. It should rarely be called directly. Also, in essence, it is very close to `SciMLBase.ODEFunction`. - ## Fields $(TYPEDFIELDS) ## Methods - `residual(pb::BifFunction, x, p)` calls `pb.F(x,p)` -- `jacobian(pb::BifFunction, x, p)` calls `pb.J(x, p)` +- `jacobian(pb::BifFunction, x, p)` calls `pb.J(x, p)` - `dF(pb::BifFunction, x, p, dx)` calls `pb.dF(x,p,dx)` - etc """ @@ -94,6 +93,7 @@ $(TYPEDFIELDS) ## Methods +- `reMake(pb; kwargs...)` modify a bifurcation problem - `getu0(pb)` calls `pb.u0` - `getParams(pb)` calls `pb.params` - `getLens(pb)` calls `pb.lens` @@ -104,8 +104,8 @@ $(TYPEDFIELDS) - `isSymmetric(pb)` calls `isSymmetric(pb.prob)` ## Constructors - -- `BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)` and `kwargs` are the fields above. You can pass your own jacobian with `J` (see [`BifFunction`](@ref) for description of the jacobian function) and jacobian adjoint with `Jᵗ`. For example, this can be used to provide finite differences based jacobian. +- `BifurcationProblem(F, u0, params, lens)` all derivatives are computed using ForwardDiff. +- `BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)` and `kwargs` are the fields above. You can pass your own jacobian with `J` (see [`BifFunction`](@ref) for description of the jacobian function) and jacobian adjoint with `Jᵗ`. For example, this can be used to provide finite differences based jacobian using `BifurcationKit.finiteDifferences`. """ struct $op{Tvf, Tu, Tp, Tl <: Lens, Tplot, Trec} <: AbstractAllJetBifProblem diff --git a/src/Results.jl b/src/Results.jl index 20d89d43..472c1415 100644 --- a/src/Results.jl +++ b/src/Results.jl @@ -23,12 +23,13 @@ $(TYPEDFIELDS) # Associated methods - `length(br)` number of the continuation steps +- `show(br)` display information about the branch - `eigenvals(br, ind)` returns the eigenvalues for the ind-th continuation step - `eigenvec(br, ind, indev)` returns the indev-th eigenvector for the ind-th continuation step -- `br[k+1]` gives information about the k-th step. A typical run yields the following - `getNormalForm(br, ind)` compute the normal form of the ind-th points in `br.specialpoint` - `getLens(br)` return the parameter axis used for the branch - `getLenses(br)` return the parameter two axis used for the branch when 2 parameters continuation is used (Fold, Hopf, NS, PD) +- `br[k+1]` gives information about the k-th step. A typical run yields the following ``` julia> br[1] (x = 0.0, param = 0.1, itnewton = 0, itlinear = 0, ds = -0.01, θ = 0.5, n_unstable = 2, n_imag = 2, stable = false, step = 0, eigenvals = ComplexF64[0.1 - 1.0im, 0.1 + 1.0im], eigenvecs = ComplexF64[0.7071067811865475 - 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.7071067811865475im 0.0 - 0.7071067811865475im]) @@ -53,14 +54,17 @@ julia> br.param - `getParams(br)` Parameters passed to continuation and used in the equation `F(x, par) = 0`. - `setParam(br, p0)` set the parameter value `p0` according to `::Lens` for the parameters of the problem `br.prob` - `getLens(br)` get the lens used for the computation of the branch +- `continuation(br, ind)` performs automatic branch switching (aBS) from ind-th bifurcation point. Typically branching from equilibrium to equilibrium, or periodic orbit to periodic orbit. +- `continuation(br, ind, lens2)` performs two parameters `(getLens(br), lens2)` continuation of the ind-th bifurcation point. +- `continuation(br, ind, probPO::AbstractPeriodicOrbitProblem)` performs aBS from ind-th bifurcation point (which must be a Hopf bifurcation point) to branch of periodic orbits. """ @with_kw_noshow struct ContResult{Tkind <: AbstractContinuationKind, Tbr, Teigvals, Teigvec, Biftype, Tsol, Tparc, Tprob, Talg} <: AbstractResult{Tkind, Tprob} "holds the low-dimensional information about the branch. More precisely, `branch[i+1]` contains the following information `(recordFromSolution(u, param), param, itnewton, itlinear, ds, θ, n_unstable, n_imag, stable, step)` for each continuation step `i`.\n - `itnewton` number of Newton iterations - - `itlinear` total number of linear iterations during corrector + - `itlinear` total number of linear iterations during newton (corrector) - `n_unstable` number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation) - - `n_imag` number of eigenvalues with positive real part and non zero imaginary part for each continuation step (to detect Hopf bifurcation). - - `stable` stability of the computed solution for each continuation step. Hence, `stable` should match `eig[step]` which corresponds to `branch[k]` for a given `k`. + - `n_imag` number of eigenvalues with positive real part and non zero imaginary part at current continuation step (useful to detect Hopf bifurcation). + - `stable` stability of the computed solution for each continuation step. Hence, `stable` should match `eig[step]` which corresponds to `branch[k]` for a given `k`. - `step` continuation step (here equal `i`)" branch::StructArray{Tbr} @@ -70,16 +74,16 @@ julia> br.param "Vector of solutions sampled along the branch. This is set by the argument `saveSolEveryStep::Int64` (default 0) in [`ContinuationPar`](@ref)." sol::Tsol - "The parameters used for the call to `continuation` which produced this branch. Must be a ContinationPar" + "The parameters used for the call to `continuation` which produced this branch. Must be a [`ContinuationPar`](@ref)" contparams::Tparc "Type of solutions computed in this branch." kind::Tkind = EquilibriumCont() - "Structure associated to the functional, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here." + "Bifurcation problem used to compute the branch, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here." prob::Tprob = nothing - "A vector holding the set of detected bifurcation points. See [`SpecialPoint`](@ref) for a description of the fields." + "A vector holding the set of detected bifurcation points. See [`SpecialPoint`](@ref) for a list of special points." specialpoint::Vector{Biftype} "Continuation algorithm used for the computation of the branch" @@ -115,11 +119,22 @@ getfirstusertype(br::AbstractBranchResult) = keys(br.branch[1])[1] setParam(br::AbstractBranchResult, p0) = setParam(br.prob, p0) Base.getindex(br::ContResult, k::Int) = (br.branch[k]..., eigenvals = haseigenvalues(br) ? br.eig[k].eigenvals : nothing, eigenvecs = haseigenvector(br) ? br.eig[k].eigenvecs : nothing) Base.lastindex(br::ContResult) = length(br) + +""" +$(SIGNATURES) + +Return the solution for the ind-th point stored in br.sol +""" @inline function getSolx(br::ContResult, ind::Int) @assert hassolution(br) "You did not record the solution in the branch. Please set `saveSolEveryStep` in `ContinuationPar`" return br.sol[ind].x end +""" +$(SIGNATURES) + +Return the parameter for the ind-th point stored in br.sol +""" @inline function getSolp(br::ContResult, ind::Int) @assert hassolution(br) "You did not record the solution in the branch. Please set `saveSolEveryStep` in `ContinuationPar`" return br.sol[ind].p @@ -143,9 +158,9 @@ Return the eigenvalues of the ind-th continuation step. `verbose` is used to tel function eigenvals(br::AbstractBranchResult, ind::Int, verbose::Bool = false) @assert br.eig[ind+1].step == ind "Error in indexing eigenvalues. Please open an issue on the website." if verbose - println("--> For ", getLensSymbol(br), " = ", br.branch[ind].param) - println("--> There are ", br.branch[ind].n_unstable, " unstable eigenvalues") - println("--> Eigenvalues for continuation step ", br.eig[ind+1].step) + println("──> For ", getLensSymbol(br), " = ", br.branch[ind].param) + println("──> There are ", br.branch[ind].n_unstable, " unstable eigenvalues") + println("──> Eigenvalues for continuation step ", br.eig[ind+1].step) end ~br.eig[ind+1].converged && @error "Eigen solver did not converged on the step!!" br.eig[ind+1].eigenvals diff --git a/src/codim2/NormalForms.jl b/src/codim2/NormalForms.jl index 49add332..1b11904f 100644 --- a/src/codim2/NormalForms.jl +++ b/src/codim2/NormalForms.jl @@ -545,15 +545,15 @@ function bogdanovTakensNormalForm(_prob, # Al-Hdaibat, B., W. Govaerts, Yu. A. Kuznetsov, and H. G. E. Meijer. “Initialization of Homoclinic Solutions near Bogdanov--Takens Points: Lindstedt--Poincaré Compared with Regular Perturbation Method.” SIAM Journal on Applied Dynamical Systems 15, no. 2 (January 2016): 952–80. https://doi.org/10.1137/15M1017491. ########################### - vext = real.(ζs[1]) + vr = real.(ζs[1]) Lᵗ = hasAdjoint(prob_vf) ? jad(prob_vf, x0, parbif) : transpose(L) _λ★, _ev★, _ = eigsolver(Lᵗ, nev) Ivp = sortperm(_λ★, by = abs) # in case the prob is HopfMA, we real it zerov = real.(prob_ma.zero) - wext = real.(geteigenvector(eigsolver, _ev★, Ivp[1])) - q0, = bls(L, wext, vext, zero(𝒯), zerov, one(𝒯)) - p1, = bls(Lᵗ, vext, wext, zero(𝒯), zerov, one(𝒯)) + vl = real.(geteigenvector(eigsolver, _ev★, Ivp[1])) + q0, = bls(L, vl, vr, zero(𝒯), zerov, one(𝒯)) + p1, = bls(Lᵗ, vr, vl, zero(𝒯), zerov, one(𝒯)) q1, = bls(L, p1, q0, zero(𝒯), q0, zero(𝒯)) p0, = bls(Lᵗ, q0, p1, zero(𝒯), p1, zero(𝒯)) # we want diff --git a/src/periodicorbit/NormalForms.jl b/src/periodicorbit/NormalForms.jl index 452b91c6..7d961521 100644 --- a/src/periodicorbit/NormalForms.jl +++ b/src/periodicorbit/NormalForms.jl @@ -116,7 +116,7 @@ function perioddoublingNormalForm(pbwrap::WrapPOSh, Teigvec = vectortype(br), detailed = true, kwargs_nf...) - verbose && println("━"^53*"\n──> Period-doubling normal form computation") + verbose && println("━"^53*"\n──▶ Period-doubling normal form computation") bifpt = br.specialpoint[ind_bif] bptype = bifpt.type pars = setParam(br, bifpt.param) @@ -249,7 +249,7 @@ function perioddoublingNormalForm(pbwrap::WrapPOColl, prm = true, kwargs_nf...) # first, get the bifurcation point parameters - verbose && println("━"^53*"\n──> Period-Doubling normal form computation") + verbose && println("━"^53*"\n──▶ Period-Doubling normal form computation") bifpt = br.specialpoint[ind_bif] bptype = bifpt.type par = setParam(br, bifpt.param) @@ -342,7 +342,7 @@ function neimarksackerNormalForm(pbwrap::WrapPOColl, coll = pbwrap.prob N, m, Ntst = size(coll) @assert coll isa PeriodicOrbitOCollProblem "Something is wrong. Please open an issue on the website" - verbose && println("━"^53*"\n──> Period-doubling normal form computation") + verbose && println("━"^53*"\n──▶ Period-doubling normal form computation") # bifurcation point bifpt = br.specialpoint[ind_bif] @@ -373,16 +373,16 @@ function neimarksackerNormalForm(pbwrap::WrapPOColl, J[end, 1:end-1] .= q J[1:end-1, end] .= p - wext = J' \ rhs - vext = J \ rhs + vl = J' \ rhs + vr = J \ rhs - v₁ = @view vext[1:end-1] - v₁★ = @view wext[1:end-1] + v₁ = @view vr[1:end-1] + v₁★ = @view vl[1:end-1] - ζₛ = getTimeSlices(coll, vext) - vext ./= sqrt(∫(coll, ζₛ, ζₛ, 1)) + ζₛ = getTimeSlices(coll, vr) + vr ./= sqrt(∫(coll, ζₛ, ζₛ, 1)) - ζ★ₛ = getTimeSlices(coll, wext) + ζ★ₛ = getTimeSlices(coll, vl) v₁★ ./= 2∫(coll, ζ★ₛ, ζₛ, 1) v₁ₛ = getTimeSlices(coll, vcat(v₁,0)) @@ -483,7 +483,7 @@ function neimarksackerNormalForm(pbwrap::WrapPOSh{ <: ShootingProblem }, # first, get the bifurcation point parameters sh = pbwrap.prob @assert sh isa ShootingProblem "Something is wrong. Please open an issue on the website" - verbose && println("━"^53*"\n──> Neimark-Sacker normal form computation") + verbose && println("━"^53*"\n──▶ Neimark-Sacker normal form computation") # bifurcation point bifpt = br.specialpoint[ind_bif] diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index c43aa513..213d7e9c 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/src/periodicorbit/PeriodicOrbitCollocation.jl @@ -665,6 +665,7 @@ end # this function updates the section during the continuation run @views function updateSection!(prob::PeriodicOrbitOCollProblem, x, par; stride = 0) + @debug "Update section Collocation" # update the reference point prob.xπ .= 0 diff --git a/src/periodicorbit/PeriodicOrbitTrapeze.jl b/src/periodicorbit/PeriodicOrbitTrapeze.jl index eaa12b2d..522a9e68 100644 --- a/src/periodicorbit/PeriodicOrbitTrapeze.jl +++ b/src/periodicorbit/PeriodicOrbitTrapeze.jl @@ -648,6 +648,7 @@ end # this function updates the section during the continuation run @views function updateSection!(prob::PeriodicOrbitTrapProblem, x, par; stride = 0) + @debug "Update section TRAP" M, N = size(prob) xc = getTimeSlices(prob, x) T = extractPeriodFDTrap(prob, x) diff --git a/src/periodicorbit/StandardShooting.jl b/src/periodicorbit/StandardShooting.jl index 00423f9e..fab2de09 100644 --- a/src/periodicorbit/StandardShooting.jl +++ b/src/periodicorbit/StandardShooting.jl @@ -99,6 +99,7 @@ end # this function updates the section during the continuation run function updateSection!(sh::ShootingProblem, x, par) + @debug "Update section shooting" xt = getTimeSlices(sh, x) @views update!(sh.section, vf(sh.flow, xt[:, 1], par), xt[:, 1]) sh.section.normal ./= norm(sh.section.normal) diff --git a/test/test_wave.jl b/test/test_wave.jl index 56e778dd..7f21aa0a 100644 --- a/test/test_wave.jl +++ b/test/test_wave.jl @@ -67,11 +67,11 @@ Fcgl(u, p, t = 0) = Fcgl!(similar(u), u, p) end #################################################################################################### n = 50 - l = pi +l = pi - Δ, D = Laplacian1D(n, l, :Periodic) - par_cgl = (r = 0.0, μ = 0.5, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ), Db = blockdiag(D, D), γ = 0.0, δ = 1.0, N = 2n) - sol0 = zeros(par_cgl.N) +Δ, D = Laplacian1D(n, l, :Periodic) +par_cgl = (r = 0.0, μ = 0.5, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ), Db = blockdiag(D, D), γ = 0.0, δ = 1.0, N = 2n) +sol0 = zeros(par_cgl.N) # _sol0 = zeros(2n) # _J0 = Jcgl(_sol0, par_cgl)