Skip to content

Commit

Permalink
add matrix based NS/PD continuation
Browse files Browse the repository at this point in the history
  • Loading branch information
@rveltz
rveltz committed Oct 7, 2024
1 parent a6e81fa commit 5be3e9a
Show file tree
Hide file tree
Showing 3 changed files with 213 additions and 92 deletions.
110 changes: 80 additions & 30 deletions src/periodicorbit/codim2/MinAugNS.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,37 +73,8 @@ end
res = 𝐍𝐒(x[begin:end-2], x[end-1], x[end], params)
return vcat(res[1], res[2], res[3])
end

###################################################################################################
# Struct to invert the jacobian of the pd MA problem.
struct NSLinearSolverMinAug <: AbstractLinearSolver; end

function NSMALinearSolver(x, p::𝒯, ω::𝒯, 𝐍𝐒::NeimarkSackerProblemMinimallyAugmented, par,
duu, dup, duω;
debugArray = nothing) where 𝒯
################################################################################################
# debugArray is used as a temp to be filled with values used for debugging.
# If debugArray = nothing, then no debugging mode is entered.
# If it is AbstractArray, then it is populated
################################################################################################
# Recall that the functional we want to solve is [F(x,p), σ(x,p)]
# where σ(x,p) is computed in the above functions and F is the periodic orbit
# functional. We recall that N⋅[v, σ] ≡ [0, 1]
# The Jacobian Jpd of the functional is expressed at (x, p)
# We solve here Jpd⋅res = rhs := [rhsu, rhsp, rhsω]
# The Jacobian expression of the NS problem is
# ┌ ┐
# Jhopf = │ J dpF 0 │
# │ σx σp σω │
# └ ┘
########## Resolution of the bordered linear system ########
# J * dX + dpF * dp = du => dX = x1 - dp * x2
# The second equation
# <σx, dX> + σp * dp + σω * dω = du[end-1:end]
# thus becomes
# (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
# This 2 x 2 system is then solved to get (dp, dω)
########################## Extraction of function names ########################################
function _get_bordered_terms(𝐍𝐒::NeimarkSackerProblemMinimallyAugmented, x, p::𝒯, ω::𝒯, par) where 𝒯
a = 𝐍𝐒.a
b = 𝐍𝐒.b

Expand Down Expand Up @@ -144,6 +115,78 @@ function NSMALinearSolver(x, p::𝒯, ω::𝒯, 𝐍𝐒::NeimarkSackerProblemMi
σω = -(dot(w, apply(jacobian_neimark_sacker(POWrap, x, par, ω+ϵ2), v)) -
dot(w, apply(jacobian_neimark_sacker(POWrap, x, par, ω), v)) )/ϵ2

return (;JNS, JNS★, dₚF, σₚ, δ, ϵ2, ϵ3, v, w, par0, dJvdp, itv, itw, σω)
end
###################################################################################################
function jacobian(pdpb::NSMAProblem{Tprob, MinAugMatrixBased}, X, par) where {Tprob}
𝐍𝐒 = pdpb.prob

# get the PO functional, ie a WrapPOSh, WrapPOTrap, WrapPOColl
POWrap = 𝐍𝐒.prob_vf

x = @view X[begin:end-2]
p = X[end-1]
ω = X[end]

𝒯 = eltype(p)

@unpack JNS★, dₚF, σₚ, ϵ2, ϵ3, v, w, par0, σω = _get_bordered_terms(𝐍𝐒, x, p, ω, par)

cw = conj(w)
vr = real(v); vi = imag(v)
u1r = apply_jacobian_neimark_sacker(POWrap, x .+ ϵ2 .* vcat(vr,0), par0, ω, cw, true)
u1i = apply_jacobian_neimark_sacker(POWrap, x .+ ϵ2 .* vcat(vi,0), par0, ω, cw, true)
u2 = apply(JNS★, cw)
σxv2r = @. -(u1r - u2) / ϵ2 # careful, this is a complex vector
σxv2i = @. -(u1i - u2) / ϵ2
σx = @. σxv2r + Complex{𝒯}(0, 1) * σxv2i

dJvdt = minus(apply(jacobian_neimark_sacker(POWrap, x .+ ϵ2 .* vcat(0 * vr, 1), par0, ω), v),
apply(jacobian_neimark_sacker(POWrap, x .- ϵ2 .* vcat(0 * vr, 1), par0, ω), v));
rmul!(dJvdt, 𝒯(1/(2ϵ3)))
σt = -dot(w, dJvdt)

_Jpo = jacobian(POWrap, x, par0)
Jns = hcat(_Jpo.jacpb, dₚF, zero(dₚF))
Jns = vcat(Jns, vcat(real(σx), real(σt), real(σₚ), real(σω))')
Jns = vcat(Jns, vcat(imag(σx), imag(σt), imag(σₚ), imag(σω))')
end
###################################################################################################
# Struct to invert the jacobian of the ns MA problem.
struct NSLinearSolverMinAug <: AbstractLinearSolver; end

function NSMALinearSolver(x, p::𝒯, ω::𝒯, 𝐍𝐒::NeimarkSackerProblemMinimallyAugmented, par,
duu, dup, duω;
debugArray = nothing) where 𝒯
################################################################################################
# debugArray is used as a temp to be filled with values used for debugging.
# If debugArray = nothing, then no debugging mode is entered.
# If it is AbstractArray, then it is populated
################################################################################################
# Recall that the functional we want to solve is [F(x,p), σ(x,p)]
# where σ(x,p) is computed in the above functions and F is the periodic orbit
# functional. We recall that N⋅[v, σ] ≡ [0, 1]
# The Jacobian Jpd of the functional is expressed at (x, p)
# We solve here Jpd⋅res = rhs := [rhsu, rhsp, rhsω]
# The Jacobian expression of the NS problem is
# ┌ ┐
# Jhopf = │ J dpF 0 │
# │ σx σp σω │
# └ ┘
########## Resolution of the bordered linear system ########
# J * dX + dpF * dp = du => dX = x1 - dp * x2
# The second equation
# <σx, dX> + σp * dp + σω * dω = du[end-1:end]
# thus becomes
# (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
# This 2 x 2 system is then solved to get (dp, dω)
########################## Extraction of function names ########################################

# get the PO functional, ie a WrapPOSh, WrapPOTrap, WrapPOColl
POWrap = 𝐍𝐒.prob_vf

@unpack JNS★, dₚF, σₚ, ϵ2, ϵ3, v, w, par0, σω, itv, itw = _get_bordered_terms(𝐍𝐒, x, p, ω, par)

if has_hessian(𝐍𝐒) == false || 𝐍𝐒.usehessian == false
cw = conj(w)
vr = real(v); vi = imag(v)
Expand Down Expand Up @@ -280,6 +323,11 @@ function continuation_ns(prob, alg::AbstractContinuationAlgorithm,
nspointguess = vcat(nspointguess.u, nspointguess.p...)
prob_ns = NSMAProblem(𝐍𝐒, FiniteDifferencesMF(), nspointguess, par, lens2, plot_solution, prob.recordFromSolution)
opt_ns_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
elseif jacobian_ma == :MinAugMatrixBased
nspointguess = vcat(nspointguess.u, nspointguess.p...)
prob_ns = NSMAProblem(𝐍𝐒, MinAugMatrixBased(), nspointguess, par, lens2, plot_solution, prob.recordFromSolution)
opt_ns_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver

else
prob_ns = NSMAProblem(𝐍𝐒, nothing, nspointguess, par, lens2, plot_solution, prob.recordFromSolution)
opt_ns_cont = @set options_cont.newton_options.linsolver = NSLinearSolverMinAug()
Expand Down Expand Up @@ -425,6 +473,8 @@ function continuation_ns(prob, alg::AbstractContinuationAlgorithm,

# eigen solver
eigsolver = HopfEig(getsolver(opt_ns_cont.newton_options.eigsolver), prob_ns)

prob_ns = re_make(prob_ns, record_from_solution = _recordsol2)

# change the plotter
_kwargs = (record_from_solution = record_from_solution(prob), plot_solution = plot_solution)
Expand Down
110 changes: 75 additions & 35 deletions src/periodicorbit/codim2/MinAugPD.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,69 @@ end
res = 𝐏𝐝(x[begin:end-1], x[end], params)
return vcat(res[1], res[2])
end
###################################################################################################
function _get_bordered_terms(𝐏𝐝::PeriodDoublingProblemMinimallyAugmented, x, p::𝒯, par) where 𝒯
a = 𝐏𝐝.a
b = 𝐏𝐝.b

# get the PO functional, ie a WrapPOSh, WrapPOTrap, WrapPOColl
POWrap = 𝐏𝐝.prob_vf

# parameter axis
lens = getlens(𝐏𝐝)
# update parameter
par0 = set(par, lens, p)

# we define the following jacobian. It is used at least 3 times below. This avoids doing 3 times the (possibly) costly building of J(x, p)
JPD = jacobian_period_doubling(POWrap, x, par0) # jacobian with period doubling boundary condition

# we do the following in order to avoid computing the jacobian twice in case 𝐏𝐝.Jadjoint is not provided
JPD★ = has_adjoint(𝐏𝐝) ? jacobian_adjoint_period_doubling(POWrap, x, par0) : transpose(JPD)

# we solve N[v, σ1] = [0, 1]
v, σ1, cv, itv = pdtest(JPD, a, b, zero(𝒯), 𝐏𝐝.zero, one(𝒯); lsbd = 𝐏𝐝.linbdsolver)
~cv && @debug "Linear solver for N did not converge."

# # we solve Nᵗ[w, σ2] = [0, 1]
w, σ2, cv, itw = pdtest(JPD★, b, a, zero(𝒯), 𝐏𝐝.zero, one(𝒯); lsbd = 𝐏𝐝.linbdsolverAdjoint)
~cv && @debug "Linear solver for Nᵗ did not converge."

δ = getdelta(POWrap)
ϵₚ = ϵₓ = ϵⱼ = ϵₜ = 𝒯(δ)

dₚF = minus(residual(POWrap, x, set(par, lens, p + ϵₚ)),
residual(POWrap, x, set(par, lens, p - ϵₚ))); rmul!(dₚF, 𝒯(1 / (2ϵₚ)))
dJvdp = minus(apply(jacobian_period_doubling(POWrap, x, set(par, lens, p + ϵⱼ)), v),
apply(jacobian_period_doubling(POWrap, x, set(par, lens, p - ϵⱼ)), v));
rmul!(dJvdp, 𝒯(1/(2ϵⱼ)))
σₚ = -dot(w, dJvdp)

return (;JPD, JPD★, dₚF, σₚ, δ, ϵₜ, ϵₓ, v, w, par0, dJvdp, itv, itw)
end
###################################################################################################
function jacobian(pdpb::PDMAProblem{Tprob, MinAugMatrixBased}, X, par) where {Tprob}
p = X[end]
x = @view X[begin:end-1]

𝐏𝐝 = pdpb.prob
𝒯 = eltype(p)

POWrap = 𝐏𝐝.prob_vf

@unpack dₚF, σₚ, ϵₜ, ϵₓ, v, w, par0 = _get_bordered_terms(𝐏𝐝, x, p, par)

u1 = apply_jacobian_period_doubling(POWrap, x .+ ϵₓ .* vcat(v,0), par0, w, true)
u2 = apply_jacobian_period_doubling(POWrap, x .- ϵₓ .* vcat(v,0), par0, w, true)
σₓ = minus(u2, u1); rmul!(σₓ, 1 / (2ϵₓ))

# a bit of a hack
xtmp = copy(x); xtmp[end] += ϵₜ
σₜ = (𝐏𝐝(xtmp, p, par0)[end] - 𝐏𝐝(x, p, par0)[end]) / (ϵₜ)

_Jpo = jacobian(POWrap, x, par0)

return [_Jpo.jacpb dₚF ; vcat(σₓ, σₜ)' σₚ]
end
###################################################################################################
# Struct to invert the jacobian of the pd MA problem.
struct PDLinearSolverMinAug <: AbstractLinearSolver; end
Expand Down Expand Up @@ -98,41 +160,11 @@ function PDMALinearSolver(x, p::𝒯, 𝐏𝐝::PeriodDoublingProblemMinimallyAu
# σx = -< w, d2F(x,p)[v, x2]>
# where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = n
########################## Extraction of function names ########################################
a = 𝐏𝐝.a
b = 𝐏𝐝.b

# get the PO functional, ie a WrapPOSh, WrapPOTrap, WrapPOColl
POWrap = 𝐏𝐝.prob_vf

# parameter axis
lens = getlens(𝐏𝐝)
# update parameter
par0 = set(par, lens, p)

# we define the following jacobian. It is used at least 3 times below. This avoids doing 3 times the (possibly) costly building of J(x, p)
JPD = jacobian_period_doubling(POWrap, x, par0) # jacobian with period doubling boundary condition

# we do the following in order to avoid computing the jacobian twice in case 𝐏𝐝.Jadjoint is not provided
JPD★ = has_adjoint(𝐏𝐝) ? jacobian_adjoint_period_doubling(POWrap, x, par0) : transpose(JPD)

# we solve N[v, σ1] = [0, 1]
v, σ1, cv, itv = pdtest(JPD, a, b, zero(𝒯), 𝐏𝐝.zero, one(𝒯); lsbd = 𝐏𝐝.linbdsolver)
~cv && @debug "Linear solver for N did not converge."

# # we solve Nᵗ[w, σ2] = [0, 1]
w, σ2, cv, itw = pdtest(JPD★, b, a, zero(𝒯), 𝐏𝐝.zero, one(𝒯); lsbd = 𝐏𝐝.linbdsolverAdjoint)
~cv && @debug "Linear solver for Nᵗ did not converge."

δ = getdelta(POWrap)
ϵₚ = ϵₓ = ϵⱼ = ϵₜ = 𝒯(δ)
################### computation of σx σp ####################
################### and inversion of Jpd ####################
dₚF = minus(residual(POWrap, x, set(par, lens, p + ϵₚ)),
residual(POWrap, x, set(par, lens, p - ϵₚ))); rmul!(dₚF, 𝒯(1 / (2ϵₚ)))
dJvdp = minus(apply(jacobian_period_doubling(POWrap, x, set(par, lens, p + ϵⱼ)), v),
apply(jacobian_period_doubling(POWrap, x, set(par, lens, p - ϵⱼ)), v));
rmul!(dJvdp, 𝒯(1/(2ϵⱼ)))
σₚ = -dot(w, dJvdp)
@unpack dₚF, σₚ, ϵₜ, ϵₓ, v, w, par0, itv, itw = _get_bordered_terms(𝐏𝐝, x, p, par)

if has_hessian(𝐏𝐝) == false || 𝐏𝐝.usehessian == false
# We invert the jacobian of the PD problem when the Hessian of x -> F(x, p) is not known analytically.
Expand Down Expand Up @@ -286,23 +318,31 @@ function continuation_pd(prob, alg::AbstractContinuationAlgorithm,
linbdsolve_adjoint = bdlinsolver_adjoint,
usehessian = usehessian)

# this is to remove this part from the arguments passed to continuation
_kwargs = (record_from_solution = record_from_solution, plot_solution = plot_solution)
_plotsol = modify_po_plot(prob, _kwargs)

@assert jacobian_ma in (:autodiff, :finiteDifferences, :minaug, :finiteDifferencesMF, :MinAugMatrixBased)

# Jacobian for the PD problem
if jacobian_ma == :autodiff
pdpointguess = vcat(pdpointguess.u, pdpointguess.p)
prob_pd = PDMAProblem(𝐏𝐝, AutoDiff(), pdpointguess, par, lens2, plot_solution, prob.recordFromSolution)
prob_pd = PDMAProblem(𝐏𝐝, AutoDiff(), pdpointguess, par, lens2, _plotsol, prob.recordFromSolution)
opt_pd_cont = @set options_cont.newton_options.linsolver = DefaultLS()
elseif jacobian_ma == :finiteDifferences
pdpointguess = vcat(pdpointguess.u, pdpointguess.p...)
prob_pd = PDMAProblem(𝐏𝐝, FiniteDifferences(), pdpointguess, par, lens2, plot_solution, prob.recordFromSolution)
prob_pd = PDMAProblem(𝐏𝐝, FiniteDifferences(), pdpointguess, par, lens2, _plotsol, prob.recordFromSolution)
opt_pd_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
elseif jacobian_ma == :finiteDifferencesMF
pdpointguess = vcat(pdpointguess.u, pdpointguess.p)
prob_pd = PDMAProblem(𝐏𝐝, FiniteDifferencesMF(), pdpointguess, par, lens2, plot_solution, prob.recordFromSolution)
prob_pd = PDMAProblem(𝐏𝐝, FiniteDifferencesMF(), pdpointguess, par, lens2, _plotsol, prob.recordFromSolution)
opt_pd_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
elseif jacobian_ma == :MinAugMatrixBased
pdpointguess = vcat(pdpointguess.u, pdpointguess.p)
prob_pd = PDMAProblem(𝐏𝐝, MinAugMatrixBased(), pdpointguess, par, lens2, _plotsol, prob.recordFromSolution)
opt_pd_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
else
prob_pd = PDMAProblem(𝐏𝐝, nothing, pdpointguess, par, lens2, plot_solution, prob.recordFromSolution)
prob_pd = PDMAProblem(𝐏𝐝, nothing, pdpointguess, par, lens2, _plotsol, prob.recordFromSolution)
opt_pd_cont = @set options_cont.newton_options.linsolver = PDLinearSolverMinAug()
end

Expand Down
Loading

0 comments on commit 5be3e9a

Please sign in to comment.