Library

Parameters

BifurcationKit.NewtonParType
struct NewtonPar{T, L<:BifurcationKit.AbstractLinearSolver, E<:AbstractEigenSolver}

Returns a variable containing parameters to affect the newton algorithm when solving F(x) = 0.

Arguments for line search (Armijo)

  • linesearch = false: use line search algorithm (i.e. Newton with Armijo's rule)
  • α = 1.0: initial value of α (damping) parameter for line search algorithm
  • αmin = 0.001: minimal value of the damping alpha
Mutating

For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package Accessors.jl to drastically simplify the mutation of different fields. See the tutorials for examples.

Internal fields

  • tol::Any: absolute tolerance for F(x). Default: 1.0e-12

  • max_iterations::Int64: number of Newton iterations. Default: 25

  • verbose::Bool: display Newton iterations? Default: false

  • linsolver::BifurcationKit.AbstractLinearSolver: linear solver, must be <: AbstractLinearSolver. Default: DefaultLS()

  • eigsolver::AbstractEigenSolver: eigen solver, must be <: AbstractEigenSolver. Default: DefaultEig()

  • linesearch::Bool: Default: false

  • α::Any: Default: convert(typeof(tol), 1.0)

  • αmin::Any: Default: convert(typeof(tol), 0.001)

source
BifurcationKit.ContinuationParType
struct ContinuationPar{T, S<:BifurcationKit.AbstractLinearSolver, E<:AbstractEigenSolver}

Returns a variable containing the parameters to affect the continuation algorithm used to solve F(x, p) = 0.

Arguments

  • dsmin, dsmax are the minimum, maximum allowed arc-length value. It controls the density of points in the computed branch of solutions.
  • ds = 0.01 is the initial arc-length.
  • p_min, p_max allowed parameter range for p
  • max_steps = 100 maximum number of continuation steps
  • newton_options::NewtonPar: options for the Newton algorithm
  • save_to_file = false: save to file. A name is automatically generated or can be defined in continuation. This requires using JLD2.
  • save_sol_every_step::Int64 = 1 at which continuation steps do we save the current solution
  • plot_every_step = 10 at which continuation steps do we plot the current solution

Handling eigen-elements, their computation is triggered by the argument detect_bifurcation (see below)

  • nev = 3 number of eigenvalues to be computed. It is automatically increased to have at least nev unstable eigenvalues. To be set for proper bifurcation detection. See Detection of bifurcation points of Equilibria for more information.
  • save_eig_every_step = 1 record eigen vectors every specified steps. Important for memory limited resource, e.g. GPU.
  • save_eigenvectors = true Important for memory limited resource, e.g. GPU.

Handling bifurcation detection

  • tol_stability = 1e-10 lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
  • detect_fold = true detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)
  • detect_bifurcation::Int ∈ {0, 1, 2, 3} If set to 0, nothing is done. If set to 1, the eigen-elements are computed. If set to 2, the bifurcations points are detected during the continuation run, but not located precisely. If set to 3, a bisection algorithm is used to locate the bifurcations points (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory of memory limited devices (e.g. on GPU)
  • dsmin_bisection = 1e-16 minimal ds for the bisection algorithm for locating bifurcation points
  • n_inversion = 2 number of sign inversions in bisection algorithm
  • max_bisection_steps = 25 maximum number of bisection steps
  • tol_bisection_eigenvalue = 1e-16 tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps

Handling ds adaptation (see continuation for more information)

  • a = 0.5 aggressiveness factor. It is used to adapt ds in order to have a number of newton iterations per continuation step roughly constant. The higher a is, the larger the step size ds is changed at each continuation step.

Handling event detection

  • detect_event::Int ∈ {0, 1, 2} If set to 0, nothing is done. If set to 1, the event locations are sought during the continuation run, but not located precisely. If set to 2, a bisection algorithm is used to locate the event (slower).
  • tol_param_bisection_event = 1e-16 tolerance on parameter to locate event

Misc

  • η = 150. parameter to estimate tangent at first point with parameter p₀ + ds / η
  • detect_loop [WORK IN PROGRESS] detect loops in the branch and stop the continuation
Mutating

For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package Accessors.jl to drastically simplify the mutation of different fields. See tutorials for more examples.

source

Problems

GridapBifurcationKit.GridapBifProblemType
GridapBifProblem(res, u0, parms, V, U, dΩ, lens; jac = nothing, mass = nothing, kwargs...)

Construct a bifurcation problem which encodes a system of PDEs discretized with Gridap. It is a subtype of BifurcationKit.AbstractDAEBifProblem so that a (possibly singular) mass matrix can be used for the stability analysis, e.g. for Hopf bifurcations.

Arguments

  • res(u, p, v): residual of the (semi-)discretized problem, where p are the parameters.
  • u0: initial guess (a Gridap FEFunction or its free dof values).
  • parms: the set of parameters.
  • V: TestFESpace.
  • U: TrialFESpace.
  • : the Measure used to assemble the residual (stored for the mass matrix).
  • lens: an Accessors lens selecting the continuation parameter in parms, e.g. (@optic _.λ).

Keyword arguments

  • jac(u, p, du, v): analytical jacobian. If nothing, it is computed by finite differences.
  • d2res(u, p, du1, du2, v), d3res(u, p, du1, du2, du3, v): second and third derivatives, required for automatic branch switching with a non-simple kernel.
  • mass(u, v): integrand of the mass matrix, e.g. (u,v) -> ∫(u⋅v)*dΩ. If nothing, the default L² mass ∫(u⋅v)*dΩ is used (see get_mass_matrix).
  • record_from_solution, plot_solution, R01, R02, R11, delta: see the BifurcationKit documentation.

Extended methods

  • get_mass_matrix assembles the mass matrix associated to mass.
  • BifurcationKit.is_mass_matrix_constant and BifurcationKit.getmassmatrix are specialized so that the mass matrix is used by the DAE eigensolvers.
source
GridapBifurcationKit.get_mass_matrixFunction
get_mass_matrix(prob::GridapBifProblem)
get_mass_matrix(prob::GridapProblem, dΩ = prob.dΩ)

Assemble the (sparse) mass matrix associated to the problem, on the free dofs of the trial/test spaces (Dirichlet dofs are eliminated, consistently with the jacobian).

The integrand is the one passed to GridapBifProblem(...; mass = ...), defaulting to the L² mass ∫(u⋅v)*dΩ. For an incompressible flow, a typical choice is the velocity only mass (u,p),(v,q) -> ∫(v⊙u)*dΩ, which yields a singular mass matrix with a zero pressure block, as required for the stability of the differential-algebraic system $M\dot z = F(z, p)$.

source

Eigen solvers

The generic eigensolvers are provided by BifurcationKit.jl (DefaultEig, EigArpack, EigArnoldiMethod, EigKrylovKit) and the DAE wrappers by EigenDAE and EigenMassMatrix. See Eigen Solvers.

BifurcationKit.EigenMassMatrixType
struct EigenMassMatrix{Tb, Teig<:AbstractEigenSolver} <: AbstractEigenSolver

Create an eigensolver for DAE, Basically a GEV with mass matrix.

Internal fields

  • B::Any: Mass matrix

  • eig::AbstractEigenSolver: Eigen-solver

source

Branch switching (branch point)

Automatic branch switching at a branch point is performed with

continuation(br::ContResult, ind_bif::Int, optionsCont::ContinuationPar; kwargs...)

where br is a branch computed with detection of bifurcation points enabled. See Branch switching for more information and the precise method definition in BifurcationKit.jl.

Branch switching (Hopf point)

Automatic branch switching at a Hopf point towards periodic orbits is performed with

continuation(br, ind_bif::Int, _contParams::ContinuationPar,	prob::AbstractPeriodicOrbitProblem ; δp = nothing, ampfactor = 1, kwargs...)

See Branch switching for more information and the precise method definition in BifurcationKit.jl.

Normal form

The normal forms of simple branch points and Hopf points are computed with the BifurcationKit function get_normal_form, see Simple bifurcation branch point and Simple Hopf point.