allow to do codim 2 continuation of PO with collocation and mesh adap…

…tation

src/periodicorbit/codim2/MinAugNS.jl

@@ -11,6 +11,20 @@ function ns_point(br::AbstractBranchResult, index::Int)
     return BorderedArray(_copy(specialpoint.x), [specialpoint.param, ω])
 end
 
+function ns_point(br::AbstractResult{Tkind, Tprob}, index::Int) where {Tkind <: PeriodicOrbitCont, Tprob <: WrapPOColl}
+    bptype = br.specialpoint[index].type
+    @assert bptype == :ns "This should be a NS point"
+    specialpoint = br.specialpoint[index]
+    ω = imag(br.eig[specialpoint.idx].eigenvals[specialpoint.ind_ev])
+    if specialpoint.x isa NamedTuple
+        # the solution is mesh adapted, we need to restore the mesh.
+        pbwrap = deepcopy(br.prob)
+        update_mesh!(pbwrap.prob, specialpoint.x._mesh )
+        specialpoint = @set specialpoint.x = specialpoint.x.sol
+    end
+    return BorderedArray(_copy(specialpoint.x), [specialpoint.param, ω])
+end
+
 function apply_jacobian_neimark_sacker(pb, x, par, ω, dx, _transpose = false)
     if _transpose == false
         @assert 1==0

src/periodicorbit/codim2/PeriodicOrbitCollocation.jl

@@ -56,6 +56,32 @@ function continuation(br::AbstractResult{Tkind, Tprob},
     nothing
 end
 
+function foldpoint(br::AbstractResult{Tkind, Tprob}, index::Int) where {Tkind <: PeriodicOrbitCont, Tprob <: WrapPOColl}
+    bptype = br.specialpoint[index].type
+    @assert bptype == :bp || bptype == :nd || bptype == :fold "This should be a Fold / BP point"
+    specialpoint = br.specialpoint[index]
+    if specialpoint.x isa NamedTuple
+        # the solution is mesh adapted, we need to restore the mesh.
+        pbwrap = deepcopy(br.prob)
+        update_mesh!(pbwrap.prob, specialpoint.x._mesh )
+        specialpoint = @set specialpoint.x = specialpoint.x.sol
+    end
+    return BorderedArray(_copy(specialpoint.x), specialpoint.param)
+end
+
+function pd_point(br::AbstractResult{Tkind, Tprob}, index::Int) where {Tkind <: PeriodicOrbitCont, Tprob <: WrapPOColl}
+    bptype = br.specialpoint[index].type
+    @assert bptype == :pd "This should be a PD point"
+    specialpoint = br.specialpoint[index]
+    if specialpoint.x isa NamedTuple
+        # the solution is mesh adapted, we need to restore the mesh.
+        pbwrap = deepcopy(br.prob)
+        update_mesh!(pbwrap.prob, specialpoint.x._mesh )
+        specialpoint = @set specialpoint.x = specialpoint.x.sol
+    end
+    return BorderedArray(_copy(specialpoint.x), specialpoint.param)
+end
+
 """
 $(SIGNATURES)
 
@@ -77,6 +103,14 @@ function continuation_coll_fold(br::AbstractResult{Tkind, Tprob},
     biftype = br.specialpoint[ind_bif].type
     bifpt = br.specialpoint[ind_bif]
 
+    # if mesh adaptation, we need to extract the solution specifically
+    if bifpt.x isa NamedTuple
+        # the solution is mesh adapted, we need to restore the mesh.
+        pbwrap = deepcopy(br.prob)
+        update_mesh!(pbwrap.prob, bifpt.x._mesh )
+        bifpt = @set bifpt.x = bifpt.x.sol
+    end
+
     # we get the collocation problem
     coll = getprob(br).prob
 
@@ -142,7 +176,7 @@ function continuation_coll_pd(br::AbstractResult{Tkind, Tprob},
 
     # get the PD eigenvectors
     par = setparam(br, bifpt.param)
-    jac = jacobian(br.prob, bifpt.x, par)
+    jac = jacobian(br.prob, pdpointguess.u, par)
     J = jac.jacpb
     nj = size(J, 1)
     J[end, :] .= rand(nj) # must be close to kernel
@@ -189,9 +223,7 @@ function continuation_coll_ns(br::AbstractResult{Tkind, Tprob},
                     kwargs...) where {Tkind <: PeriodicOrbitCont, Tprob <: WrapPOColl}
     bifpt = br.specialpoint[ind_bif]
     biftype = bifpt.type
-
     @assert biftype == :ns "We continue only NS points of Periodic orbits for now"
-
     nspointguess = ns_point(br, ind_bif)
 
     # we copy the problem for not mutating the one passed by the user
@@ -200,10 +232,10 @@ function continuation_coll_ns(br::AbstractResult{Tkind, Tprob},
 
     # get the NS eigenvectors
     par = setparam(br, bifpt.param)
-    jac = jacobian(br.prob, bifpt.x, par)
+    jac = jacobian(br.prob, nspointguess.u, par)
     J = Complex.(copy(jac.jacpb))
     nj = size(J, 1)
-    J[end, :] .= rand(nj) # must be close to eigensapce
+    J[end, :] .= rand(nj) # must be close to eigenspace
     J[:, end] .= rand(nj)
     J[end, end] = 0
     # enforce NS boundary condition

src/periodicorbit/codim2/utils.jl

@@ -1,63 +0,0 @@
-# function modify_po_2params_finalise(prob, kwargs, probMA)
-#     update_section_every_step = prob.update_section_every_step
-#     _finsol = get(kwargs, :finalise_solution, nothing)
-#     if isnothing(_finsol)
-#         return (Z, tau, step, contResult; kF...) ->
-#             begin
-#                 z = Z.u
-#                 # we first check that the continuation step was successful
-#                 # if not, we do not update the problem with bad information
-#                 success = converged(get(kF, :state, nothing)) && ~get(kF, :bisection, false)
-#                 if success && mod_counter(step, update_section_every_step) == 1
-#                     updatesection!(prob, getvec(z, probMA), setparam(contResult, Z.p))
-#                 end
-#                 return true
-#             end
-#     else
-#         return (Z, tau, step, contResult; kF...) ->
-#             begin
-#                 # we first check that the continuation step was successful
-#                 # if not, we do not update the problem with bad information!
-#                 success = converged(get(kF, :state, nothing)) && ~get(kF, :bisection, false)
-#                 if success && mod_counter(step, update_section_every_step) == 1
-#                     z = Z.u
-#                     updatesection!(prob, getvec(z, probMA), setparam(contResult, Z.p))
-#                 end
-#                 return _finsol(Z, tau, step, contResult; prob = prob, kF...)
-#             end
-#     end
-# end
-
-# function modify_po_2params_finalise(prob::PeriodicOrbitOCollProblem, kwargs, probMA)
-#     update_section_every_step = prob.update_section_every_step
-#     _finsol = get(kwargs, :finalise_solution, nothing)
-#     _finsol2 = (Z, tau, step, contResult; kF...) ->
-#         begin
-#             # we first check that the continuation step was successful
-#             # if not, we do not update the problem with bad information
-#             success = converged(get(kF, :state, nothing)) && ~get(kF, :bisection, false)
-#             # mesh adaptation
-#             z = Z.u
-#             if success && prob.meshadapt
-#                 oldsol = _copy(z)
-#                 oldmesh = getTimes(prob) .* getPeriod(prob, getvec(z, probMA), nothing)
-#                 adapt = computeError!(prob, getvec(z, probMA);
-#                         verbosity = prob.verboseMeshAdapt,
-#                         par = setparam(contResult, z.p),
-#                         K = prob.K)
-#                 if ~adapt.success
-#                     return false
-#                 end
-#             end
-#             if success && mod_counter(step, update_section_every_step) == 1
-#                 # parameter is not used for section in collocation
-#                 updatesection!(prob, getvec(z, probMA), nothing)
-#             end
-#             if isnothing(_finsol)
-#                 return true
-#             else
-#                 return _finsol(Z, tau, step, contResult; prob = prob, kF...)
-#             end
-#         end
-#     return _finsol2
-# end