External solver interface

Residuals and derivatives for solver integration.

See the API overview for other analyses and conventions for units and array ordering.

JosephsonCircuits.HBNonlinearProblem Type

HBNonlinearProblem

The harmonic balance system as a nonlinear problem: the residual, its exact derivatives and the pieces a preconditioner reads, with the unknowns a plain real vector so that a solver written elsewhere can drive it.

Build one with hbnonlinearproblem. length(p) is the number of unknowns; hbresidual!, hbjvp!, hbjacobian! and hbvjp! evaluate at any point, and setdrive! scales the drive for continuation.

When the circuit injects direct current the unknowns are the canonical state rather than the harmonic one, and augmentation carries the block which makes the difference. isaugmented says which, and every entry point above routes through the block when it is there, so a caller does not have to know. See DCAugmentation.

Fields

• sys: the harmonic balance system the residual is evaluated on.

• modelayout: the layout relating the real unknowns to the complex modes.

• u0: the point the problem was built at.

• jacobian: the assembled Jacobian, or nothing when none was asked for.

• parts: the pieces a preconditioner needs.

• bnm0: the drive as built, which setdrive! scales.

• tplan: the transposed gather maps the vector-Jacobian product walks.

• Pwork, Qwork, betawork, dirtd3: transform workspaces.

• augmentation: the direct current block, or nothing.

• atol: the residual tolerance a root is held to, norm(F) <= atol with the norm finite: the solve's atol, raised to the rounding floor of the source, or its rtol times the residual at u0 when that is finite, whichever is larger.

source

JosephsonCircuits.JacobianOperator Type

JacobianOperator(prob::HBNonlinearProblem, u)

The Jacobian of prob at u as a linear operator.

Implements size, eltype, LinearAlgebra.mul!, and the same for its adjoint and transpose (which are matrix free through the transposed gather maps; see hbvjp!). Construction sets the system's evaluation point ONCE – one forward transform – and neither mul! touches it again unless something else has moved it since, so the products inside a Krylov loop, in either direction, pay exactly two transforms each and nothing more. The operator keeps a copy of u, so moving u afterwards, in place or not, leaves it where it was; construct a new operator at the new point: that is the point update, and it is the same cost the internal solver pays once per Newton step.

J = JacobianOperator(prob, u)
Krylov.gmres(J, -F; N = preconditioner(prob, u), atol = 0.0)

Set the absolute tolerance to zero inside a Newton loop

rtol in a Krylov solver is relative to the norm of the right hand side, and in a Newton loop that right hand side is the residual being driven to zero. Any absolute floor eventually exceeds it, at which point the linear solver correctly reports success for a system it never touched, the Newton step is zero, and the iteration stagnates with nothing reporting a failure.

Krylov.jl defaults to atol = sqrt(eps()), about 1.5e-8, which is sensible standalone and wrong here. This package's own gmres! defaults to atol = 0.0, and nlsolvekrylov! passes a tenth of its nonlinear tolerance, tying the floor to that rather than to machine epsilon.

source

JosephsonCircuits.drivenresidual! Method

drivenresidual!(F, prob, u, scale)

The residual at u with the drive scaled by scale: the F(u, p) form a continuation library asks for, with p the drive scale.

source

JosephsonCircuits.hbd2F! Method

hbd2F!(out, prob, u, v, w)

The exact second directional derivative H(u)[v,w] = B(f''(A*u).*(A*v).*(A*w)), with f'' the second derivative of the current-phase relation of each junction, -sin for a Josephson junction; the linear term contributes nothing.

source

JosephsonCircuits.hbd3F! Method

hbd3F!(out, prob, u, v, w, z)

The exact third directional derivative d3F(u)[v,w,z] = B(f'''(A*u).*(A*v).*(A*w).*(A*z)), with f''' the third derivative of the current-phase relation of each junction, -cos for a Josephson junction.

Supplying this exactly rather than by finite differences is what makes a continuation library's normal form computation – cusp, Bogdanov-Takens, Bautin – accurate.

source

JosephsonCircuits.hbdFdp! Method

hbdFdp!(out, prob)

The derivative of the residual with respect to the drive scale, which is -b: the residual is B(f(A*u)) + K*u - scale*b, with f the current-phase relation of each junction, so the parameter derivative is exact and constant, and a continuation tangent never needs a finite difference in the parameter.

source

JosephsonCircuits.hbjacobian! Method

hbjacobian!(J, prob::HBNonlinearProblem, u)

Assemble the exact real Jacobian at u.

source

JosephsonCircuits.hbjvp! Method

hbjvp!(Jv, prob::HBNonlinearProblem, u, v)

The exact matrix-free Jacobian-vector product J(u)*v. Two transforms plus the linear term, with no Jacobian assembled.

This convenience form re-sets the evaluation point on every call, which costs one extra forward transform. Inside a Krylov loop use a JacobianOperator, which sets the point once at construction and whose mul! pays only the product itself.

source

JosephsonCircuits.hbnonlinearproblem Function

hbnonlinearproblem(w, Nharmonics, sources, circuit,
    circuitdefs = Dict{Symbol,Any}(); assemblejacobian = true, kwargs...)

Build the harmonic balance system as an HBNonlinearProblem without solving it, by calling hbnlsolve with returnsystem = true and the remaining keywords; circuitdefs is optional when every component value is numeric, as it is there. The problem lives on the host: a device backend is refused. assemblejacobian = false leaves out the assembled real Jacobian, which a matrix-free solver does not need. The tolerance a solve would hold a root to, from atol and rtol, is the problem's atol.

source

JosephsonCircuits.hbresidual! Method

hbresidual!(F, prob::HBNonlinearProblem, u)

The harmonic balance residual at u, in place.

source

JosephsonCircuits.hbvjp! Method

hbvjp!(out, prob, u, w)

The transposed product transpose(J(u))*w, matrix free.

Needed by BifurcationKit for minimally augmented codim-2 continuation, fold and Hopf continuation and left eigenvectors, and accepted by NonlinearSolve as vjp.

The adjoints of the two transforms compose so that the normalizations cancel, leaving the conjugate multiplicity divided out in one kernel and multiplied back in through the coefficients of the other; see NonlinearTermTransposePlan.

source

JosephsonCircuits.preconditioner Method

preconditioner(prob::HBNonlinearProblem, u; spec = BlockDiagonal(),
    precision = Float64)

The mode coupling preconditioner of prob (see ModeCouplingPreconditioner for spec, a member of the mode coupling family), updated at u, its matrix assembled in precision, in which it is factorized too, except that KLU and UMFPACK factorize in double precision whatever they are handed and a factorization with a precision of its own holds its factors in that (factorizationprecision). With an explicit direct current block it is wrapped in the canonical coordinates with the direct current subsystem solved exactly.

Applied by ldiv! and by mul!, so it can be handed straight to any external Krylov solver: Krylov.gmres(J, -F; N = preconditioner(prob, u), atol = 0.0). See JacobianOperator for why atol = 0.0. updatepreconditioner!(P, u) refactorizes it at another point, reusing its structure and symbolic analysis. Building or updating it sets the problem's evaluation point, which an operator built at another point notices and sets back.

source

JosephsonCircuits.setdrive! Method

setdrive!(prob::HBNonlinearProblem, scale)

Scale the drive of prob by scale, in place, and return the problem.

The residual is F(u) = B(f(A*u)) + K*u - b, with f the current-phase relation of each junction (sin for a Josephson junction, the polynomial of a PolynomialCPR for a nonlinear inductor), and the drive enters only through b. Scaling it is therefore the one parameter which can be varied without rebuilding anything: the sparsity, the plans and the preconditioner structure are all untouched.

This is what makes continuation in pump power possible. scale = 0 is the undriven problem, whose solution is zero and whose Jacobian is the linear circuit; stepping up to scale = 1 and carrying the converged state forward walks onto the driven branch, which is the reliable way to reach an operating point a cold solve cannot find.

The scale is relative to the drive the problem was built with, so setdrive!(prob, 1) restores it.

source