Harmonic-balance analyses

Solve operating points and signal responses, reuse setup, and differentiate outputs.

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

JosephsonCircuits.ContourIntegral Type

julia

ContourIntegral(center, radius; quadrature = 64, moments = 4,
    probes = 12, ranktol = 1e-10, changetol = 1e-6, refinements = 3,
    tol = 1e-9)

The method of hbstability which finds the poles inside the circle of center and radius, in inverse seconds, by the block moment integral of Beyn (Linear Algebra Appl. 436, 2012, 3839), the method for LaplaceResponse models and exact TransmissionLine delays, whose operator is not a polynomial; Monodromy takes the lines of a pumped circuit as well. It factors the operator at quadrature points of the circle, refactoring its fixed pattern, and doubles the points, keeping those it has, up to refinements times, until the poles move by less than changetol of the radius and the moments by less than changetol of the mean norm of the solves sampled on the circle. probes right hand sides and moments moments bound how many poles it holds, ranktol cuts the moments' rank relative to the larger of their largest singular value and that mean norm, so that an empty circle converges, and a pole is accepted with a residual under tol. A saturated rank or a failed check sets converged = false; the checks are numerical, not a certified count of the poles. Keep the poles away from the circle, and the branch cuts of an analytic response outside it.

JosephsonCircuits.DenseSpectrum Type

julia

DenseSpectrum(; tol = 1e-9, maxunknowns = 1000)

The method of hbstability which finds every finite pole of the truncated problem, its Floquet aliases included, by the QZ algorithm on its companion pencil formed densely, twice the unknowns square: the reference for a small circuit. Its time grows as the cube of the unknowns, every node's harmonics counted, so a problem of more than maxunknowns is refused before its pencil is formed; Monodromy finds the selected least damped modes of a pumped circuit, and ShiftInvert the poles near given frequencies. A pole is accepted with a residual under tol.

JosephsonCircuits.HB Type

julia

HB(nonlinear, linearized)

The result of hbsolve: the NonlinearHB solution for the pump and its harmonics in nonlinear, and the LinearizedHB solution for the signals in linearized.

JosephsonCircuits.HBCache Type

julia

HBCache

A reusable harmonic balance solver over a typed circuit with parameters: the compiled circuit and a copy of its definitions, the mode grid and its Fourier index maps, the solver options, and the last converged operating point, which hbsolve! uses to warm start the next solve.

Built by hbcache. converged reports whether the last point converged, from its warm start or from the cold start a failed warm start is retried from (see hbsolve!), and a point which does not converge also warns with the reason it stopped and leaves the stored operating point as it was. Check it: a solve which does not converge returns a state that looks like a solution and is not one, and comparing timings or gradients against it is meaningless. A cache is mutable and every solve through it rewrites what it holds, so it belongs to one solve at a time: do not share a cache between concurrent solves; build one per task.

JosephsonCircuits.HBStabilityResult Type

julia

HBStabilityResult

The finite poles hbstability found, with time dependence exp(s*t), in decreasing order of their real parts, every array in that order:

poles: the growth rates s in inverse seconds, growing where the real part is positive.
residuals: each pole's relative backward error in the scaled polynomial, which meets the method's tol, or under Monodromy with a pump a condition-based estimate of multiplier roundoff relative to the multiplier. It excludes timestep error, which rateerrors estimates separately.
rateerrors: for pumped Monodromy, an estimate of each real part's timestep error. The selected modes are propagated through a period in twice the steps, their line histories resampled and read back in the coarse coordinates, and the refined map projected onto their right and left vectors, G = (VL' VR) \ VL' M2 VR, which follows their mixing among themselves. Each mode continues into the eigenvector of G that carries most of it, of multiplier μ2; the estimate is (16/15*abs(log(abs(μ2/μ))) + residuals[k])/T, the coarse map's error where it falls as the fourth power of the step. Inf means no refined mode continues the mode. Mixing with modes beyond nev is not followed, and a map too coarse for its error to fall so can have more error than its estimate. abs(real(poles[k])) > rateerrors[k] is evidence of a resolved sign, growth where real(poles[k]) > rateerrors[k] and decay where real(poles[k]) < -rateerrors[k], not a certificate: confirm it by doubling steps, which cuts the estimate about sixteenfold once the steps resolve the mode, and where modes crowd by a larger nev. nothing for other methods and without a pump.
edgeweights: the fraction of each mode's squared amplitude, node voltage over frequencyscale and junction flux together, in the two outermost harmonics, zero without a pump. A large one says to retain more harmonics in a polynomial method, and a wider profile under the period map; a small one bounds no truncation error.
nodevoltage, junctionflux: each mode's coefficients, (mode, node, pole) in volts and (mode, junction, pole) in webers, normalized together, along modes, nodes (ground excluded) and junctionbranches (the junctions' branches in the compiled circuit, whose orientation sets the sign). For a nonzero profile, sum(abs2, nodevoltage[:,:,k]/frequencyscale) + sum(abs2, junctionflux[:,:,k]) is one. The mode's amplitude and global phase are arbitrary; these coefficients do not predict a driven oscillation amplitude. Harmonic n has angular frequency imag(s) + n*pumpfrequency.
stateedgeweights: the same fraction for the states of rational blocks, which depends on their realization, and under Monodromy the waves leaving the lines' ports with them, in the voltage they carry over frequencyscale.
internalonly: modes whose node/junction content is negligible beside their balanced rational-state content or, under pumped Monodromy, the waves leaving their lines' ports; their exported profiles are zero and their dominant harmonic is selected from those states and waves. Zero profiles do not imply that a mode is absent or stable.
harmonics: each mode's dominant harmonic, the n whose coefficients carry most of its squared amplitude, by the norm of edgeweights (by its states' and waves' for a mode internal to blocks and lines): the mode lives mostly at the angular frequency imag(s) + n*pumpfrequency, which tells a pole's Floquet aliases apart by their content, and one whose dominant harmonic is an outermost retained one is shaped by the truncation.
shiftindices: the shift of ShiftInvert each pole was found near, zero for the other methods.
searches: diagnostics for each search. A shift records requested and converged counts and products taken. A contour records quadrature, rank and capacity, saturation, changes in poles and moments, and the nearest pole's distance to the circle. Pumped Monodromy records steps, period, the multipliers not classified as unresolved (multipliers), how many of the largest were tested against their own condition-based bounds (tested), gauge directions removed, and line-history columns. Multipliers within eps()*opnorm(M,1) of zero after balancing are classified as unresolved without an individual test and do not count in tested. With nev = :all every multiplier above that threshold is tested against its own bound; with a smaller nev the untested ones are taken against the threshold alone, as well conditioned, so the counts can depend on nev. The array is empty for DenseSpectrum and the unpumped Monodromy fallback.
converged: whether the finite eigensolve/search met its acceptance checks. Polynomial/contour candidates must meet tol; rejected counts failed candidates. Pumped Monodromy classifies unresolved multipliers separately and can return converged=true with Inf rate errors. This field does not establish discretization accuracy, search completeness for the physical circuit, or stability.
infinite: the infinite eigenvalues DenseSpectrum set aside, those ShiftInvert's searches came upon (see there), or the multipliers of the period map within the eigensolve's error of zero, eps times the norm of the balanced map, and those of the largest within their own, that over their condition: its algebraic directions and the modes damped beyond a period's precision. The map tests its largest multipliers against their own bounds until nev pass, and the others against the bound of a well conditioned one alone, so that the count is the eigensolve's for every multiplier with nev = :all and can be smaller with fewer; nothing for ContourIntegral.
pumpfrequency, frequencyscale: angular frequencies in rad/s.
method: the method which found the poles.

A physical mode has Floquet aliases s + im*k*pumpfrequency. Harmonic truncation can perturb these copies differently; polynomial methods return the candidates they find. A pole found near two shifts is returned twice. The period map returns one representative per selected mode, up to nev. converged concerns the eigensolve: neither it nor a dense spectrum establishes that the circuit has no unstable mode the truncation in harmonics or in steps misses.

JosephsonCircuits.LinearizedHB Type

julia

LinearizedHB(w, modes, S, Snoise, Cnoise, Vout, Ssensitivity, QE,
    QEideal, CM, nbar, nodeflux, nodefluxadjoint, voltage,
    voltageadjoint, nodenames, nodeindices, componentnames,
    componenttypes, componentnamedict, mutualinductorbranchnames,
    portnumbers, portindices, portimpedances, porttemperatures,
    noiseportimpedanceindices, channeltemperatures, sensitivitynames,
    sensitivityindices, Nmodes, Nnodes, Nbranches, Nports, signalindex)

The solution of the linearized harmonic balance problem returned by hblinsolve. An output which was not requested is an empty array. As plain arrays the scattering parameters and the quantum efficiencies are indexed [output, input, frequency] and the sensitivities [output, input, component, frequency], with the mode and the port of each flattened (mode fastest); with keyedarrays = true they are keyed arrays with the axes outputmode, outputport, inputmode, inputport, component and freqindex. CM, nbar, Snoise, Cnoise, Vout and the node outputs have axes of their own, named in their keys.

Fields

w: the signal angular frequencies in radians per second.
modes: the retained signal modes as tuples of harmonic indices; (0,) is the signal itself and (k,) the idler offset by k pump harmonics.
S: the scattering matrix relating the inputs and outputs at each combination of port, mode and signal frequency.
Snoise: the scattering matrix from the noise channels of the dissipative elements to the ports, in the same signed frequency convention as S: a resistor's channel is a port of its impedance in vacuum, and its row of Snoise is what that port's column of S would be (see adjointnoisesigns!). The channels are the dissipative lumped components (see noiseindices) followed by one per port of each dissipative ScatteringParameters. Being a scattering matrix it describes a transformation and does not depend on temperature.
Cnoise: the symmetrized noise covariance the circuit adds at the output ports, in quanta, sum_c (nbar_c + 1/2)*Snoise[c,i]*conj(Snoise[c,j]) over its noise channels c, when returnCnoise = true. A vacuum channel counts as 1/2, the vacuum's half photon, and a channel at temperature T as nbar + 1/2, nbar its thermaloccupation. Together with S this is the Gaussian channel the circuit implements, taking the symmetrized covariance of the inputs sigma, I/2 for the vacuum, to S*sigma*S' + Cnoise; the temperature of each channel enters here. It is the noise the circuit adds and never that of its ports, so NoiseCovariance(Cnoise) embeds the circuit as a block. Both of its port mode indices are outputs: keyed, its axes are outputmode and outputport for i, conjoutputmode and conjoutputport for j, and freqindex (see Cnoisetokeyed).
Vout: the symmetrized covariance of the waves leaving the ports, S*Diagonal(sigma)*S' + Cnoise with sigma the noise every input brings, nbar + 1/2 at the temperature of its port's termination (see MatchedTermination), when returnVout = true. Its axes are those of Cnoise, and its diagonal is nbar + 1/2.
Ssensitivity: the derivative of S with respect to a relative (logarithmic) perturbation of each component in sensitivitynames, or of each design parameter when the sensitivity pair interface is used, at a fixed pump operating point or including the shift of the operating point: by default from hbsolve (sensitivityoperatingpoint), and from hblinsolve when it is given sensitivityresidual.
QE: the quantum efficiency at each combination of port, mode and frequency: abs2(S[i,j])/(2*Vout[i,i]), the photon gain over twice the noise at the output, with every input in its state, the input j included, so a warm source lowers the efficiency of a measurement of its signal.
QEideal: the quantum efficiency of an ideal amplifier with the same gain.
CM: the bosonic commutation relations of each output, sum_j |S[i,j]|^2 sign(w_j) over the input modes, which equal +1 for an output at positive frequency and -1 for one at negative frequency when the scattering matrix is complete.
nbar: the occupation of each output port mode, the mean number of photons per mode of the wave leaving the port, Vout[i,i] - 1/2: zero in the vacuum, and the photons an amplifier emits from the vacuum it amplifies included. At an input port it is what the circuit sends back toward the device there. Its axes are those of CM.
nodeflux: the node fluxes resulting from a unit current source at each port and mode.
nodefluxadjoint: the node fluxes of the adjoint problem, the transposed linearized system driven by a unit current source at each port and mode, which the noise outputs read (see assemblesystemmatrix!).
voltage: the node voltages resulting from a unit current source at each port and mode.
voltageadjoint: the node voltages of the adjoint problem.
nodenames: the node names, ground first.
nodeindices: the 2 by Ncomponents matrix of component node indices from the CompiledCircuit.
componentnames, componenttypes, componentnamedict: the corresponding fields of the CompiledCircuit.
mutualinductorbranchnames: the names of the two inductors each mutual inductor couples, see coupledinductornames.
portnumbers: the port numbers, in the order the port axes use.
portindices: the flat component index of each port.
portimpedances: the reference impedance of each port, which the scattering parameters are normalized to.
porttemperatures: the physical temperature in kelvin of each port's termination, in the order of portnumbers.
noiseportimpedanceindices: the flat component indices of the internal dissipative components, in the order of the noise channel axis of Snoise.
channeltemperatures: the physical temperature in kelvin of each noise channel, in the order of the channel axis of Snoise: with Snoise it gives each channel's share of the noise at an output, (thermaloccupation(w, T) + 1/2)*abs2(Snoise[c, i]); a channel of a block which states its noise is at zero.
sensitivitynames: the component names the sensitivities were taken with respect to.
sensitivityindices: their flat component indices.
Nmodes: the number of retained signal modes.
Nnodes: the number of nodes, including ground.
Nbranches: the number of branches in the circuit graph.
Nports: the number of ports.
signalindex: the position of the signal mode (0,) in modes, which is always 1.

JosephsonCircuits.Monodromy Type

julia

Monodromy(; steps = 64, nev = 10, backend = CPU())

Find the least damped modes of a periodic operating point with hbstability. The method forms a dense period map, computes all of its multipliers, and returns up to nev numerically resolved modes (10 by default). nev = :all requests all resolved modes. Each selected mode has one frequency representative rather than the aliases returned by a harmonic polynomial search.

The HB orbit initializes transientsolve, which records one pump period T in steps Gauss–Legendre steps. transienttangent then propagates the constrained circuit perturbations: node fluxes and rates, rational-block states, and outgoing-wave histories of transmission lines. The map satisfies δx(T) = M*δx(0); its eigenvalues are multipliers μ = exp(s*T). Algebraic auxiliary currents are constrained, the starting perturbations meet the circuit's algebraic constraints, among them the voltages an ideal block (whose I + S is singular) ties between nodes with capacitance, or ties to the states of a rational or pumped block whose direct term is such a block, and uniform flux gauge directions are removed. This is a circuit state/history map, not a scattering matrix. Pumped blocks convert fully on, independent of their startup envelope, as they do in harmonic balance.

For each selected eigenvector, propagation and Fourier analysis of its periodic part choose a dominant-frequency representative of log(μ)/T and supply a mode profile. Nmodulationharmonics = (H,) sets the width of that profile only; it does not truncate the map eigenproblem. Require 2H + 1 <= steps. Frequencies outside the temporal sampling bandwidth can be placed at aliases; increase steps to resolve them. A line adds about 2*delay/(T/steps) history coordinates, and its delay must be at least one step. A mode with no node or junction content, internal to rational blocks or to lines, takes its representative from its block states and the waves leaving its lines' ports.

Refine both the HB orbit's harmonics and steps. rateerrors estimates each rate's timestep error from the selected modes propagated through a period of twice the steps: the refined map projected onto them, which follows their mixing among themselves, gives each mode's change of rate, scaled by 16/15 to the coarse map's error at the rule's fourth order. It is not a certified bound: mixing with modes beyond nev is not followed, and a map too coarse for its error to fall as the fourth power of the step can have more error than its estimate. Compare solves at increasing steps, and at a larger nev where modes crowd, before deciding the sign of a small rate. The Gauss–Legendre rule is not L-stable: decay much faster than a step can appear as slow discrete decay. Resolve the circuit's fast decays as well as its high frequencies. Delay-history interpolation also adds error.

backend selects propagation and the dense eigensolve:

CPU() computes all multipliers through a real Schur decomposition, then forms only the selected left/right vectors on the session's BLAS threads.
CUDABackend() requires CUDA.jl and CUDSS.jl, and a CUDA runtime with cuSOLVER 11.7.1 or newer for its geev, which CUDA.jl's wrapper checks at the eigensolve. The device computes all right vectors and obtains selected left vectors from their basis, each checked against its residual; a singular or ill-conditioned basis falls back to the CPU's Schur vectors. Balancing, selection, orbit initialization, and constraint setup stay on the host.

For map dimension d, storage is O(d^2) and the dense eigensolve is O(d^3) even for small nev. Reducing nev saves selected-vector and profile work. More steps also enlarge the map when it contains delays.

Without a pump this method uses DenseSpectrum() instead, with its default 1000-unknown limit. Exact unpumped delays and any LaplaceResponse model require ContourIntegral. A frequency-domain block without a supported transient realization requires an appropriate polynomial or contour method.

JosephsonCircuits.NonlinearHB Type

julia

NonlinearHB(w, frequencies, nodeflux, Rbnm, Ljb, Lb, Ljbm, Nmodes,
    Nbranches, nodes, ports, modes, S, solverinfo, operatingpoint,
    dcnodevoltage, sources, blockcurrents)

The solution of the nonlinear harmonic balance problem returned by hbnlsolve.

Fields

w: the tuple of pump angular frequencies in radians per second, one per non-commensurate pump.
frequencies: the Frequencies describing the retained pump harmonics and intermodulation products.
nodeflux: the node fluxes at each retained mode. With keyedarrays = true an Nmodes by Nnodes - 1 keyed array with axes :outputmode and :node; otherwise a vector of length Nmodes*(Nnodes - 1) with the mode index varying fastest.
Rbnm: the incidence matrix between the branch and node bases, with each entry repeated Nmodes times.
Ljb: sparse vector of the Josephson junction inductances by branch.
Lb: sparse vector of the linear inductances by branch.
Ljbm: Ljb with each entry repeated Nmodes times.
Nmodes: the number of retained modes.
Nbranches: the number of branches in the circuit graph.
nodes: the node names, ground ("0") first.
ports: the port numbers.
modes: the retained modes as tuples of harmonic indices; (1,) is the pump of a single pump solve, (1,0) the first of two pumps.
S: the scattering matrix at the pump frequencies, relating the inputs and outputs at each combination of port and mode. S[i,j] is the output wave of the solution at the port mode i divided by the incident wave at the driven port mode j, so a column is the column of the scattering matrix only when its port mode is the only one a source drives; with sources at several port modes, as commensurate pumps are given (see hbsolve), each column holds the response to all of them, and the columns of undriven port modes are zero. Its zero frequency entries are identically zero: the waves are in units of sqrt(photons/second), whose normalization 1/sqrt(|w|) (see portwavescale) has no limit at zero, so there is no direct current wave to report. The direct current operating point is in dcnodevoltage.
solverinfo: diagnostics of the solution process; see SolverInfo.
operatingpoint: the converged operating point with the exact real Jacobian there, when the solver was called with returnoperatingpoint = true; otherwise nothing. See HBOperatingPoint.
dcnodevoltage: the average voltage of each non-ground node in volts, keyed by node name like nodeflux, when the analysis has a zero frequency mode; otherwise nothing. This is distinct from the zero mode of nodeflux, which is the static flux setting the inductor currents and junction phases. A vector of zeros is an answer (a node shorted to ground sits at zero volts, and so does every node of a circuit into which no direct current is injected) and differs from nothing.
sources: the drive of the solution, the sources as named tuples (mode, port, current) with a complex current: at a nonzero mode the current in time is 2 real(current*cis(f*t)), f the mode's frequency, and at the zero mode the current itself.
blockcurrents: the port currents of the scattering blocks in amperes, entering each block at its signal terminals, as the solve determined them with the whole circuit: one matrix per block of the compiled circuit, in its order (a TransmissionLine is one), indexed [mode, port] with the modes of modes; in time they are the series nodeflux is.

JosephsonCircuits.ShiftInvert Type

julia

ShiftInvert(shifts; nev = 6, krylovdim = max(30, 3nev), restarts = 200,
    tol = 1e-9)

The method of hbstability which finds the nev poles nearest each of the complex shifts, in inverse seconds; imag(shift) = 2pi*f searches near the frequency f in Hz. Each shift takes one factorization of the polynomial there, and ArnoldiMethod's restarted Schur iteration on the companion pencil shifted and inverted through it, which never forms the pencil, with at most krylovdim basis vectors and restarts restarts, to the tolerance tol, which bounds the residual a pole is accepted with as well. The iteration starts from a fixed vector, and where its Krylov space closes before it holds the poles asked for, as it does at a repeated pole, continues from vectors drawn from a generator seeded alike at every search, so a search gives the same poles at every call and leaves the caller's random stream alone. A shift exactly at a pole makes the polynomial singular there; move it off by a little.

A search for more poles than lie near its shift comes upon the pencil's infinite eigenvalues, which the inversion maps to zero. Of its candidates which fail the polynomial, those which outnumber the finite poles it can still have are infinite: the degree of the polynomial's determinant, which the sums of its rows' and of its columns' degrees bound, less the poles the search accepted. The others are rejected, a remote finite pole the search did not resolve among them, and the search has not converged.

JosephsonCircuits.designjacobian Method

julia

designjacobian(circuit, values, p::NamedTuple; jacobian = nothing)

The Jacobian of the component values of a typed Circuit, or of its compiled circuit, with respect to the design parameters p, a NamedTuple of real numbers. values(p) returns the values of the names the circuit's components are written in as a NamedTuple of numbers, the definitions a solve of the circuit takes; a derived value is computed there, in Julia. The derivatives of the values are jacobian(p) when it is given, a matrix with a row per name of values(p) and a column per parameter of p, in their orders, and otherwise those of ForwardDiff through values, which is loaded for them (using ForwardDiff), whose real and imaginary parts are differentiated as real functions of the real parameters, so a value a parameter rotates in the complex plane is exact. A frequency dependent value is no value of values(p): no design parameter may move it, and it is written into the circuit as the component's value.

Returns (names, values, J): the component names in compiled order, their values, and the complex length(names) by length(p) Jacobian, dense, materialized from the sparse enumeration of the nonzero derivatives (designdirections). A port's slot holds its reference impedance: a port which owns a termination is not among the components, its reference impedance being the value of that termination, which is; a port which owns none (termination = nothing) is, and a parameter moves only the normalization of its waves through it. A component whose value is frequency dependent is not among them either. A scattering block has no scalar value to differentiate; its dependence is carried separately by designblockjacobian.

JosephsonCircuits.designsensitivities Method

julia

designsensitivities(circuit, values, p::NamedTuple, ws, wp, sources,
    Nmodulationharmonics, Npumpharmonics; jacobian = nothing,
    keyedarrays = true, kwargs...)

The derivative of the scattering parameters with respect to the design parameters p, a NamedTuple of real numbers, of a typed Circuit whose component values are written as names, which the function values gives, values(p) a NamedTuple of numbers, by the chain rule through the component sensitivities:

julia

dS/dp_j = sum_k (dS/dv_k) (dv_k/dp_j).

The derivatives of the values are those of designjacobian, by ForwardDiff through values, loaded for them (using ForwardDiff), or from jacobian(p), so a derived value like C/4 contributes its factor of one quarter without being declared; a ScatteringParameters block depends on a parameter of p through the analytic derivative its derivatives states for that parameter's name (designblockjacobian). A port's reference impedance normalizes its waves, so a parameter it depends on moves that normalization too, through the termination the port owns or, for a port which owns none, through the port itself. The solve runs once, at the definitions values(p), with the adjoint sensitivities of hbsolve carrying the exact direction dv_k/dp_j of each dependent component into the contraction (sensitivityoperatingpoint = true, so the shift of the pump operating point is included). A tone frequency or a source current may be written as a name values(p) gives, which the solve takes, but the derivative carries the component values alone: a drive which a design parameter moves is an ArgumentError rather than a derivative which leaves its motion out.

Carrying the direction, rather than scaling a relative derivative, is what makes the result exact for complex component values: a parameter which rotates a value in the complex plane, a loss tangent, say, has dv/dp not parallel to v, which no single relative derivative can represent. All the components a parameter touches merge into one contraction, so a design variable shared across a long line costs one contraction rather than one per cell.

Returns (out, dSdp): the full hbsolve output, and a keyed array dS/dp with the axes of the scattering sensitivity and a parameter axis, the names of p, in place of the component axis; with keyedarrays = false, the plain array of the solve's scattering sensitivity, (output, input, parameter, frequency), the port and the mode of each side flattened as in the plain S, its component axis the parameters. Additional keyword arguments are forwarded to hbsolve, as keyedarrays is.

Extended help

A gradient based optimizer wants a closure from a parameter vector to a value and a derivative, evaluated at the same point in sequence, so memoize one solve for both. For the gain in dB, G_k = 20*log10(abs(S[out, in, k])), the chain rule is dG/dp = (20/log(10))*real(conj(S)*dSdp)/abs2(S):

julia

circuit = Circuit([(:p1, 1, 0, Port(1)), (:cc, 1, 2, Capacitor(:Cc)),
    (:jj, 2, 0, JosephsonJunction(:Lj)), (:cj, 2, 0, Capacitor(1000e-15))])
values(p) = (Lj = p.Lj, Cc = p.Cc)
mutable struct Objective; lastp; lastr; end
const obj = Objective(nothing, nothing)
function solveat(pvec)
    if obj.lastp != pvec
        obj.lastr = designsensitivities(circuit, values,
            (Lj = pvec[1], Cc = pvec[2]), ws, wp, sources, (2,), (8,))
        obj.lastp = copy(pvec)
    end
    return obj.lastr
end
value(pvec) = [20*log10(abs(s))
    for s in solveat(pvec).out.linearized.S((0,),1,(0,),1,:)]
function jacobian(pvec)
    r = solveat(pvec)
    S = r.out.linearized.S((0,),1,(0,),1,:)
    # the parameter axis precedes the frequency axis of dSdp
    d = permutedims(r.dSdp((0,),1,(0,),1,:,:))
    return (20/log(10)).*real.(conj.(S).*d)./abs2.(S)
end

JosephsonCircuits.hbcache Method

julia

hbcache(w, Nharmonics, sources, circuit, circuitdefs = Dict{Symbol,Any}();
    dc = false, odd = true, even = false, maxintermodorder = Inf,
    Nevaluationharmonics = map(i -> 2i, Nharmonics),
    frequencywindow = (0, Inf), kwargs...)

A reusable nonlinear solver over a typed Circuit, or the CompiledCircuit of one, whose values are written as names (symbols), circuitdefs giving each a number. The cache keeps a copy of circuitdefs, so editing the dictionary afterwards does not move it. hbsolve! takes the parameters to move as a named tuple keyed by their names, reads the rest from the definitions, and evaluates the values at the point; a name the definitions do not hold is refused, and so is a value which is not a number at the definitions (a frequency dependent one). The current of a source and a tone frequency may be names as well, and move with their definitions, so a sweep of a pump current, of a direct current bias or of a pump frequency reuses the cache as a sweep of a component value does. A tone frequency written as a name takes the default frequencywindow, since a window would retain other modes as the frequency moves, and the frequency of no mode may change its sign along the sweep (see HBReuse).

The harmonic selection keywords match hbnlsolve; the remaining keywords are stored and forwarded to every solve as keywords of hbnlsolve on the compiled circuit, so they are the solver keywords (method, atol, rtol, iterations, backend, ...) and are validated here: a keyword the compiled circuit solve does not accept is an ArgumentError at construction rather than a failure at the first solve, as are x0 and reuse, which the cache manages itself (the warm start through warmstart, the reuse object internally), keyedarrays = true, since the state is kept as plain vectors (false is accepted as what the cache does anyway), returnsystem = true and debugJacobian = true, which return the system rather than solve it, and method = Staged(), which the cache does not support since the continuation builds its own systems at its own truncations.

Examples

julia

circuit = Circuit([(:P1, 1, 0, Port(1; Z0 = 50.0)), (:C1, 1, 2, Capacitor(:Cc)),
    (:Lj1, 2, 0, JosephsonJunction(:Lj)), (:C2, 2, 0, Capacitor(1000e-15))])
cache = hbcache((2*pi*4.75e9,), (8,),
    [(mode=(1,), port=1, current=1e-8)], circuit,
    Dict(:Lj => 1000e-12, :Cc => 100e-15))
for Lj in (900:25:1100)*1e-12
    sol = hbsolve!(cache, (Lj = Lj,))
    cache.converged || break
end
# the pump current and frequency as parameters: the current swept upward
# from the definitions, then the frequency at the last current
cache = hbcache((:wp,), (8,),
    [(mode=(1,), port=1, current=:Ip)], circuit,
    Dict(:Lj => 1000e-12, :Cc => 100e-15, :Ip => 1e-9, :wp => 2*pi*4.75e9))
for Ip in (1:10)*1e-9
    sol = hbsolve!(cache, (Ip = Ip,))
    cache.converged || break
end
for wp in 2*pi*(4.70:0.01:4.80)*1e9
    sol = hbsolve!(cache, (wp = wp,))
    cache.converged || break
end

JosephsonCircuits.hblinsolve Function

julia

hblinsolve(w, circuit, circuitdefs; Nmodulationharmonics = (0,),
    nonlinear = nothing, threewavemixing = false,
    fourwavemixing = true, maxharmonics = Nmodulationharmonics,
    maxintermodorder = Inf, nbatches = Base.Threads.nthreads(),
    returnS = true, returnSnoise = false,
    returnCnoise = false, returnVout = false, returnQE = true,
    returnCM = true, returnnbar = true,
    returnnodeflux = false, returnnodefluxadjoint = false,
    returnvoltage = false, returnvoltageadjoint = false,
    keyedarrays = true, temperature = 0.0,
    sensitivitynames::AbstractVector = String[],
    sensitivityresidual = nothing, sensitivitymode = :auto,
    returnSsensitivity = false,
    factorization = nothing, backend = CPU())

Sweep the weak signal frequencies w through the circuit linearized about the operating point nonlinear found by hbnlsolve, or through the linear circuit when nonlinear = nothing. Any number of signal and idler modes, ports and pumps is supported. The scattering parameters, noise scattering parameters, quantum efficiency, commutation relations, node fluxes and voltages, and sensitivities are computed on request.

The linearized system is solved in the same modified nodal analysis formulation as the nonlinear one (see hbnlsolve), without gauge fixing rows because a mode at (numerically) zero total frequency is rejected with an ArgumentError; estimate a direct current limit from a sequence of decreasing nonzero frequencies instead.

Arguments

w: the signal angular frequency or frequencies in radians per second, a real number or any iterable of them.
circuit: a typed Circuit or a CompiledCircuit. A Circuit is compiled with the default node ordering; for another one, pass compile(circuit; sorting = ...).
circuitdefs: a dictionary from the symbols or symbolic variables used as component values to their numerical values. Optional when every component value is numeric.

Keywords

Nmodulationharmonics = (0,): how many harmonics of each pump to retain around the signal, which sets the signal and idler modes; (0,) is the signal alone.
nonlinear = nothing: the NonlinearHB operating point to linearize about, or nothing for a linear circuit. The sweep reads the pump modulation, the Fourier coefficients of the derivative of each junction's current-phase relation, off the operating point, and every component value, the junction inductances included, from circuitdefs. An operating point solved at other values is therefore a pump held fixed while the values move, which is what the sensitivities at a fixed operating point differentiate; for the circuit's own operating point, solve it at the same values, as hbsolve does.
threewavemixing = false: retain the odd pump harmonics around the signal, through which three wave mixing couples the modes.
fourwavemixing = true: retain the even pump harmonics around the signal, through which four wave mixing couples the modes.
maxharmonics = Nmodulationharmonics: an upper bound on the absolute harmonic index retained for each pump; see truncfreqs.
maxintermodorder = Inf: keep a mode which mixes two or more pumps only when its absolute harmonic indices sum to at most this order. Every harmonic of a single pump is kept up to its cap in maxharmonics whatever the order, so with one pump it removes nothing (see truncfreqs).
nbatches = Base.Threads.nthreads(): the number of batches the signal frequencies are split into for multithreading; 1 runs single threaded.
returnS = true: return the scattering parameters of the linearized solve.
returnSnoise = false: return the noise scattering parameters.
returnCnoise = false: return the added noise covariance at the output ports; see LinearizedHB.
returnVout = false: return the covariance of the waves leaving the ports, the inputs' noise at the temperatures of their terminations included.
returnnbar = true: return the occupation of each output port mode, the photons of the wave leaving it.
returnQE = true: return the quantum efficiency.
returnCM = true: return the commutation relations.
returnnodeflux = false, returnvoltage = false: return the node fluxes and voltages of the linearized solve.
returnnodefluxadjoint = false, returnvoltageadjoint = false: return the node fluxes and voltages of the adjoint linearized solve, the transposed system driven at each port and mode.
keyedarrays = true: return the outputs as keyed arrays with named, labeled axes rather than plain arrays.
temperature = 0.0: the physical temperature in kelvin of every dissipative element which does not state its own, and so of the noise it adds. A channel at temperature T carries the symmetrized noise nbar + 1/2, nbar its thermaloccupation, which at zero temperature is the vacuum's half photon. Raising it lowers the quantum efficiency and changes Cnoise but leaves Snoise and the commutation relations alone, since those describe the transformation rather than the state. It does not warm the ports: a port's termination states its own temperature, Port(n; termination = MatchedTermination(temperature = T)), zero unless given. A component may state its own temperature (Resistor(R; temperature = T), or a ScatteringParameters with noise = ThermalEquilibrium(T)), and takes this default otherwise. A block which states its noise with a NoiseCovariance, as an amplifier given by its scattering parameters does, adds that covariance whatever the temperature.
sensitivitynames = String[]: the components to take sensitivities with respect to, named by their identifiers as symbols or strings. Supported types are C, L, R and Lj with numeric values. A ScatteringParameters block has no scalar value to perturb and cannot be named; see designsensitivities for sensitivities with respect to the parameters of a block.
sensitivitypairs, sensitivityblockpairs, nsensitivityparameters, sensitivitylabels: the design parameter interface used by designsensitivities, which names physical parameters rather than components; each pair (componentname, parameterindex, alpha) gives the relative direction alpha = (dv/dp)/v of a component value under a parameter. Not intended to be passed directly.
sensitivityresidual = nothing: the derivatives of the harmonic balance residual with respect to each component value, the columns of calcresidualsensitivity, to include the shift of the operating point in the sensitivities; hbsolve supplies these. Without them the pump operating point is held fixed.
sensitivitymode = :auto: the order in which the operating point shift is contracted into the sensitivities. :forward costs one product against the linearized system per component and signal frequency; :reverse pushes the output functionals through the transposed pump Jacobian once per output port and mode pair, so its cost does not grow with the number of components. :auto chooses :reverse when the components outnumber the output port and mode pairs eight times over, or when the forward order's stamps, a value per stored entry of the linearized system per component, would exceed FORWARDSENSITIVITYSTAMPBYTES. Both take the parameters of scattering blocks as they take components, and support any number of pumps.
returnSsensitivity = false: return dS/dr, the derivative of the scattering matrix with respect to a relative perturbation r of each named component value (p -> r*p at r = 1), computed by the adjoint method.
factorization = nothing: the factorization of the linearized system matrix at each frequency: KLUfactorization, LUfactorization, CUDSSFactorization on a device, or BlockFactorization for dense node blocks (Nmodes unknowns per node, SparseBlockFactorization). nothing takes the backend's sparse factorization for one tone, KLU on the host and cuDSS on a device, and the block factorization for two or more tones when its factors, a set per host batch, fit in half the free memory (linearizedfactorization). The block factorization pivots only within its dense blocks and stops at a singular one, which a nonsingular matrix can have: chosen by nothing, the sweep then runs again on the backend's sparse factorization, and given explicitly it throws the SingularException. On a device the choice also picks the solver of the batch: a sparse factorization is solved by cuDSS, a BlockFactorization by the batched block factorization. QRfactorization is refused: the noise, the sensitivities and the adjoint node outputs solve the transposed system on the factors of each frequency, which the sparse QR does not provide. The precision of the solutions is the factorization's: BlockFactorization(precision = Float32) refines single precision factors against the double residual to double accuracy, and BlockFactorization(precision = Float32, refine = false) solves each frequency entirely in single precision, for the cases where single precision scattering parameters are enough, with an accuracy which falls with the conditioning of the system (see BlockFactorization); the outputs are returned in double either way.
backend = CPU(): the KernelAbstractions backend the sweep is solved on. On a device the system matrices of a batch of signal frequencies, which share one sparsity pattern, are assembled by one kernel and factorized and solved as a uniform batch: with a sparse factorization by cuDSS (CUDSSFactorization), and with a BlockFactorization, the default from two tones, by the batched dense block factorization, whose batch is sized by the free memory of the device. The adjoint (transposed) solve and the noise scattering parameters run on the device too, and the sensitivities of S on host threads from the solutions the device returns. A component whose value depends on the frequency is stamped on the device from its function's values, which the host evaluates per batch, each distinct function once a mode frequency. The sweep falls back to the host when sensitivities with respect to scattering block parameters are requested, whose stamps are rebuilt per frequency on the host; a scattering block whose parameters cannot be evaluated on the device forms its noise channels on the host from the whole adjoint solution copied back.
symfreqvar = nothing: deprecated, the parameter a frequency dependent value was written as an expression in. Write the value as a FrequencyDependent closure of the frequency instead.
returnZ, returnZadjoint, returnZsensitivity, returnZsensitivityadjoint: removed; passing any of them warns. Compute impedances from the scattering parameters instead.

Returns

LinearizedHB: the linearized solution; see LinearizedHB.

Examples

julia

circuit = Circuit(
    [:p1 => Port(1; Z0 = :Rleft),
     :l1 => Inductor(:Lm),
     :l2 => Inductor(:Lm),
     :k1 => MutualInductor(:K1, :l1, :l2),
     :cc => Capacitor(:Cc),
     :jj3 => JosephsonJunction(:Lj),
     :jj4 => JosephsonJunction(:Lj),
     :cj => Capacitor(:Cj),
     :gnd => Ground()],
    [[(:p1, 1), (:l1, 1), (:cc, 1)],
     [(:cc, 2), (:l2, 1), (:jj4, 1), (:cj, 1)],
     [(:l2, 2), (:jj3, 1)],
     [(:p1, 2), (:l1, 2), (:jj3, 2), (:jj4, 2), (:cj, 2), (:gnd, 1)]])
circuitdefs = Dict{Symbol,Complex{Float64}}(
    :Lj =>2000e-12,
    :Lm =>10e-12,
    :Cc => 200.0e-15,
    :Cj => 900e-15,
    :Rleft => 50.0,
    :Rright => 50.0,
    :K1 => 0.9,
)

Idc = 1e-6*0
Ip=5.0e-6
wp=2*pi*5e9
ws=2*pi*5.2e9
# modulation settings
Npumpharmonics = (16,)
Nmodulationharmonics = (2,)
threewavemixing=false
fourwavemixing=true

nonlinear=hbnlsolve(
    (wp,),
    Npumpharmonics,
    [
        (mode=(0,),port=1,current=Idc),
        (mode=(1,),port=1,current=Ip),
    ],
    circuit,circuitdefs;dc=true,odd=fourwavemixing,even=threewavemixing)

linearized = JosephsonCircuits.hblinsolve(ws,
    circuit, circuitdefs; Nmodulationharmonics = Nmodulationharmonics,
    nonlinear = nonlinear, threewavemixing=false,
    fourwavemixing=true, returnnodeflux=true, keyedarrays = false)
isapprox(linearized.nodeflux[:, :, 1],
    ComplexF64[9.9710247e-12 - 6.4969574e-14im 4.1038803e-15 - 5.4827649e-17im 1.8800941e-14 - 2.4486235e-16im;
     4.1012274e-15 - 1.5739758e-16im 1.0122555e-11 - 1.9579258e-13im -1.3339106e-15 + 5.1109092e-17im;
     1.8801186e-14 + 2.2521165e-16im -1.3347983e-15 - 1.5588056e-17im 9.9346965e-12 + 5.9535252e-14im;
     -7.5189204e-12 + 4.9043116e-14im 9.429854e-13 - 6.4338888e-15im 5.1933438e-12 - 3.3701835e-14im;
     7.5108044e-14 - 1.441184e-15im 1.5880415e-12 - 3.071641e-14im -3.0163552e-14 + 5.7628299e-16im;
     6.7046559e-12 + 3.9870563e-14im -4.9349771e-13 - 2.786924e-15im -2.2848852e-11 - 1.3700897e-13im;
     -1.6478932e-11 + 1.0742503e-13im 9.3800228e-13 - 6.3757553e-15im 5.1701213e-12 - 3.3440622e-14im;
     7.0277033e-14 - 1.2778916e-15im -7.5160529e-12 + 1.4537689e-13im -2.8570197e-14 + 5.2285147e-16im;
     6.6811977e-12 + 3.9629044e-14im -4.9182556e-13 - 2.7702393e-15im -3.1763778e-11 - 1.90433e-13im],
    rtol = 1e-6)

# output
true

JosephsonCircuits.hblinsolve Method

julia

hblinsolve(w, psc::CompiledCircuit, circuitdefs,
    signalfreq::Frequencies; nonlinear = nothing, kwargs...)

The linearized sweep on an already compiled circuit psc, at the signal mode set signalfreq. circuitdefs is the usual dictionary; in its place the vector of resolved component values nm.vvn of a CircuitMatrices already built for the circuit is accepted, which is what hbsolve passes so that the values are resolved once; the sweep w is a real number or any iterable of them (sweepfrequencies). This is what the other methods call after building those; it takes every keyword of the general method except the ones which describe the mode set (Nmodulationharmonics, threewavemixing, fourwavemixing, maxharmonics, maxintermodorder). Among them are the design parameter sensitivity keywords described under hbsolve (sensitivitypairs, sensitivityblockpairs, nsensitivityparameters, sensitivitylabels), which the general method passes on, and debuglsys = false, which returns the HBLinearizedSystem and its ingredients instead of solving, for building reference implementations in tests.

Examples

julia

circuit = Circuit(
    [:p1 => Port(1; Z0 = :Rleft),
     :l1 => Inductor(:Lm),
     :l2 => Inductor(:Lm),
     :k1 => MutualInductor(:K1, :l1, :l2),
     :cc => Capacitor(:Cc),
     :jj3 => JosephsonJunction(:Lj),
     :jj4 => JosephsonJunction(:Lj),
     :cj => Capacitor(:Cj),
     :gnd => Ground()],
    [[(:p1, 1), (:l1, 1), (:cc, 1)],
     [(:cc, 2), (:l2, 1), (:jj4, 1), (:cj, 1)],
     [(:l2, 2), (:jj3, 1)],
     [(:p1, 2), (:l1, 2), (:jj3, 2), (:jj4, 2), (:cj, 2), (:gnd, 1)]])
circuitdefs = Dict{Symbol,Complex{Float64}}(
    :Lj =>2000e-12,
    :Lm =>10e-12,
    :Cc => 200.0e-15,
    :Cj => 900e-15,
    :Rleft => 50.0,
    :Rright => 50.0,
    :K1 => 0.9,
)

Idc = 1e-6*0
Ip = 5.0e-6
wp = 2*pi*5e9
ws = 2*pi*5.2e9
Npumpharmonics = (2,)
Nmodulationharmonics = (2,)
threewavemixing = false
fourwavemixing = true

frequencies = JosephsonCircuits.removeconjfreqs(
    JosephsonCircuits.truncfreqs(
        JosephsonCircuits.calcfreqsrdft(Npumpharmonics),
        dc = true, odd = true, even = false, maxintermodorder = Inf,
    )
)
fi = JosephsonCircuits.fourierindices(frequencies)
Nmodes = length(frequencies.modes)
psc = JosephsonCircuits.compile(circuit)
nm = JosephsonCircuits.numericmatrices(psc, circuitdefs, Nmodes = Nmodes)
nonlinear = hbnlsolve(
    (wp,),
    [
        (mode=(0,),port=1,current=Idc),
        (mode=(1,),port=1,current=Ip),
    ],
    frequencies, fi, psc, nm)
signalfreq =JosephsonCircuits.truncfreqs(
    JosephsonCircuits.calcfreqsdft(Nmodulationharmonics),
    dc = true, odd = threewavemixing, even = fourwavemixing,
    maxintermodorder = Inf,
)
linearized = JosephsonCircuits.hblinsolve(ws, psc, circuitdefs,
    signalfreq;nonlinear = nonlinear, returnnodeflux=true, keyedarrays = false)
isapprox(linearized.nodeflux[:, :, 1],
    ComplexF64[9.970981e-12 - 6.4969003e-14im 4.0916792e-15 - 5.4667846e-17im 1.8880731e-14 - 2.4587299e-16im;
     4.0890344e-15 - 1.5692383e-16im 1.0122555e-11 - 1.9579255e-13im -1.2213006e-15 + 4.6803527e-17im;
     1.8880977e-14 + 2.2613553e-16im -1.2221136e-15 - 1.4289686e-17im 9.9345856e-12 + 5.9533921e-14im;
     -7.5316458e-12 + 4.9127007e-14im 9.388239e-13 - 6.4061256e-15im 5.2169479e-12 - 3.3849327e-14im;
     7.4881532e-14 - 1.4367631e-15im 1.5880496e-12 - 3.0716563e-14im -2.8139989e-14 + 5.3766382e-16im;
     6.7342779e-12 + 4.0036396e-14im -4.5607882e-13 - 2.5817544e-15im -2.2890664e-11 - 1.3725915e-13im;
     -1.6491603e-11 + 1.0750831e-13im 9.3385621e-13 - 6.3481653e-15im 5.1936261e-12 - 3.3587035e-14im;
     7.0064897e-14 - 1.2739631e-15im -7.5160441e-12 + 1.4537671e-13im -2.6679133e-14 + 4.8869432e-16im;
     6.7107197e-12 + 3.9793888e-14im -4.5454621e-13 - 2.5664472e-15im -3.1805451e-11 - 1.9068175e-13im],
    rtol = 1e-6)

# output
true

JosephsonCircuits.hbnlsolve Method

julia

hbnlsolve(w::NTuple{N,Any}, Nharmonics::NTuple{N,Int}, sources,
    circuit, circuitdefs; iterations = 1000,
    Nevaluationharmonics = map(i -> 2i, Nharmonics),
    maxintermodorder = Inf, dc = false, odd = true, even = false,
    atol = 1e-8, rtol = 0.0, method = NewtonKrylov(), x0 = nothing,
    keyedarrays = true,
    sensitivitynames = String[], returnoperatingpoint = false,
    frequencywindow = (0, Inf), backend = CPU(), debugJacobian = false,
    returnsystem = false, assemblejacobian = true,
    warnnotconverged = true)

Solve the nonlinear harmonic balance problem of a circuit driven by any number of strong tones (pumps) at any number of ports, including direct current biases and flux pumping through a current source and a mutual inductor. hblinsolve linearizes the circuit about the operating point found here; hbsolve runs the two in sequence.

The system is solved in a modified nodal analysis (MNA) formulation in the node flux basis: mutually coupled inductor branches and the ports of scattering blocks are assigned auxiliary branch current variables with their constitutive relations kept as explicit equations, which is algebraically equivalent to the nodal formulation wherever the latter is well posed. Promoting the coupled inductors keeps the system matrix entries bounded as the coupling coefficient approaches one, where the nodal inverse inductance entries diverge as 1/(1-k^2). The system is nondimensionalized by the solver inductance scale Z0/w0 (see calcsolverscale), the geometric mean port impedance over the geometric mean nonzero drive frequency, so the residual tolerance atol is independent of the unit system and the auxiliary variables have magnitudes comparable to the node fluxes. One gauge fixing equation per floating inductive or Josephson subnetwork and zero frequency mode makes circuits which are singular at direct current in a purely nodal formulation (nodes or subnetworks with no inductive path to ground) solvable without workaround inductors; if the net direct current injected into such a subnetwork is nonzero no periodic solution exists and an ArgumentError is thrown. The reported residual norms are those of the augmented system, and the returned structure contains only the node fluxes and the original incidence matrix. Commensurate drive frequencies whose retained intermodulation products reach (numerically) zero frequency are rejected with an ArgumentError. See src/circuit/mna.jl.

Arguments

w::NTuple{N,Any}: a tuple containing the angular frequencies of the strong tones (or pumps) such as (2pi5.0e9,) for a single tone at 5 GHz and (2pi5.0e9,2pi6.0e9) for a tone at 5 GHz and a tone at 6 GHz. The frequencies should be non-commensurate. For commensurate frequencies, the lowest frequency should be provided here, and the other added to sources with a mode index equal to the ratio. A frequency is a real number, or a name, as a component value may be, which circuitdefs gives a number.
Nharmonics::NTuple{N,Int}: the largest absolute harmonic index of each tone retained as an unknown, so the modes of the returned solution. The length of the tuple must equal the number of non-commensurate tones. The nonlinearity is evaluated on the larger Nevaluationharmonics grid.
sources: the sources, as named tuples with the fields mode, port and current: mode is a tuple of one harmonic index per tone, port is the integer port number and current the (complex) current amplitude. Any iterable of them is accepted and converted to the canonical SourceTuple form by sourcetable. For example: [(mode=(1,0),port=1,current=Ip1),(mode=(0,1),port=1,current=Ip2)] specifies two pumps where the frequency of the first pump would be 1wp1 + 0wp2 and the second 0wp1+1wp2 where wp1 is the first pump frequency and wp2 is the second pump frequency. Both of the pumps are applied to port 1 with currents Ip1 and Ip2, respectively. A current may be a name, as a component value may be, which circuitdefs gives a number.
circuit: a typed Circuit or a CompiledCircuit. A Circuit is compiled with the default node ordering; for another one, pass compile(circuit; sorting = ...).
circuitdefs: a dictionary from the symbols or symbolic variables used as component values or source currents to their numerical values. Optional when every component value and current is numeric.

Keywords

iterations = 1000: the maximum number of nonlinear solver iterations before it returns unconverged.
Nevaluationharmonics = map(i -> 2i, Nharmonics): the harmonics of each tone on the grid where the nonlinearity is sampled, at least Nharmonics. The nonlinear products of the retained modes reach higher orders, and a grid with M harmonics folds order p back to p - (2M + 1) (aliasing). Products of two retained modes reach twice Nharmonics and fold outside the retained set once the grid is half again as large, the three halves rule, but the leading nonlinearity of a junction is cubic, and products of three retained modes reach three times Nharmonics, which the three halves grid folds onto the tone itself. Twice the retained set is the default: it dealiases the cubic products and leaves the fifth order ones, whose folded contributions land outside the retained set. Only the transforms grow, the unknowns are the same.
maxharmonics: deprecated and ignored with a warning; Nharmonics is the retained set and Nevaluationharmonics the sampling grid.
symfreqvar = nothing: deprecated, the parameter a frequency dependent value was written as an expression in. Write the value as a FrequencyDependent closure of the frequency instead.
frequencywindow = (0, Inf): a lower and upper bound, in the units of w, on the absolute frequency abs(dot(w, mode)) of the retained modes, a truncation by frequency beside the truncations by order; the zero frequency mode follows dc. See truncfreqs for why a floor matters on incommensurate tones.
maxintermodorder = Inf: the largest intermodulation order retained. A mode which mixes two or more tones is kept when the absolute values of its harmonic indices sum to at most this order; every harmonic of a single tone is kept up to that tone's cap in Nharmonics, whatever the order, so for one tone this removes nothing.
dc = false: retain the zero frequency mode. A CurrentSource component of the netlist is a constant current, which drives this mode; a nonzero one without the mode is an error. It drives its current through itself from its first terminal to its second: it draws the current from the node at its first terminal and delivers it to the node at its second, the opposite sense of a port source, which injects its current into the port's first (positive) terminal.
odd = true: retain the odd harmonics, which four wave mixing couples through.
even = false: retain the even harmonics, which three wave mixing couples through.
atol = 1e-8: the residual tolerance norm(F) <= atol at which the nonlinear solution is considered converged. F is scaled by Z0/w0; see calcsolverscale. A tolerance below the rounding error of the scaled source, which a circuit whose interior is far from its port impedances can reach, is raised to it: no iteration resolves a residual below the rounding of the terms which make it.
rtol = 0.0: a relative residual tolerance; the nonlinear solve is converged when norm(F) <= max(atol, rtol*norm(F0)) with F0 the initial residual. It applies to the direct and Krylov methods; a Staged method ignores it, holding its final solve to atol and its interior stages to its own interioratol.
method = NewtonKrylov(): the nonlinear solver, an AbstractHBNonlinearSolver carrying its own options: NewtonKrylov (the matrix-free real Jacobian with a preconditioner, the default), Newton (the assembled real Jacobian), QuasiNewton (the complex holomorphic Jacobian, an approximation, with Anderson acceleration), Staged (source continuation on an adaptively grown harmonic grid, see stagedhbnlsolve, the strategy for operating points the direct methods fail outright; a search which ends without the point reports where it stalled, which does not prove that no operating point exists) or ExternalSolver.
x0 = nothing: an initial value for the node fluxes, either of the node flux length or of the full augmented length including the auxiliary variables of the modified nodal analysis formulation. x0, rtol, debugJacobian, returnsystem and assemblejacobian apply to the direct and Krylov methods. A Staged method builds and warm starts its own stages: it refuses x0, debugJacobian = true and returnsystem = true, and ignores rtol and assemblejacobian.
keyedarrays = true: return nodeflux and S as keyed arrays with named axes rather than plain arrays.
sensitivitynames = String[]: the components, named by their identifiers as symbols or strings, of a sensitivity calculation which uses this solve's operating point. The pump solve takes no sensitivities itself, and refuses a name the circuit does not have.
returnoperatingpoint = false: assemble and return the exact real Jacobian at the converged solution in the operatingpoint field, for sensitivities which include the shift of the operating point.
backend = CPU(): the KernelAbstractions backend the solve runs on.
debugJacobian = false: instead of solving, return a named tuple with the residual and Jacobian functions and the ingredients they are assembled from, for building reference implementations in tests.
returnsystem = false: instead of solving, return a named tuple with the HBSystem, the initial real state and residual, the real representation layout and (when assemblejacobian = true) the assembled real Jacobian, for driving an external solver.
assemblejacobian = true: assemble the real Jacobian for returnsystem and for an ExternalSolver; false skips the assembly when the external solver is matrix free.
ftol: deprecated, read as atol with a warning.
switchofflinesearchtol, alphamin: deprecated and ignored with a warning.
factorization: deprecated; beside the default method it solves with Newton(factorization = factorization) with a warning, and beside another method it is refused: the factorization is an option of the method.

A solve which does not converge returns the last iterate with solverinfo.converged = false and warns with the reason it stopped, which is also the reason of its IterationInfo (:iterations, :work, :linesearch, :progress, :nonfinite or :external; stallmessage spells each out). Check solverinfo.converged before using the result.

warnnotconverged = true: warn when the nonlinear solve does not converge, and run the checks of a converged point which warn (a junction carrying nearly its critical current at direct current). A continuation whose stage solves are expected to fail passes false, and reports and checks its own outcome.

Returns

NonlinearHB: A simple structure to hold the harmonic balance solutions. See NonlinearHB.

Examples

julia

circuit = Circuit(
    [:p1 => Port(1; Z0 = :Rleft),
     :l1 => Inductor(:Lm),
     :l2 => Inductor(:Lm),
     :k1 => MutualInductor(:K1, :l1, :l2),
     :cc => Capacitor(:Cc),
     :jj3 => JosephsonJunction(:Lj),
     :jj4 => JosephsonJunction(:Lj),
     :cj => Capacitor(:Cj),
     :gnd => Ground()],
    [[(:p1, 1), (:l1, 1), (:cc, 1)],
     [(:cc, 2), (:l2, 1), (:jj4, 1), (:cj, 1)],
     [(:l2, 2), (:jj3, 1)],
     [(:p1, 2), (:l1, 2), (:jj3, 2), (:jj4, 2), (:cj, 2), (:gnd, 1)]])
circuitdefs = Dict{Symbol,Complex{Float64}}(
    :Lj =>2000e-12,
    :Lm =>10e-12,
    :Cc => 200.0e-15,
    :Cj => 900e-15,
    :Rleft => 50.0,
    :Rright => 50.0,
    :K1 => 0.9,
)

Idc = 50e-5
Ip=0.0001e-6
wp=2*pi*5e9
Npumpmodes = 2
out=hbnlsolve(
    (wp,),
    (Npumpmodes,),
    [
        (mode=(0,),port=1,current=Idc),
        (mode=(1,),port=1,current=Ip),
    ],
    circuit,circuitdefs;dc=true,odd=true,even=false)
# each entry, the direct current fluxes and the much smaller pump ones
all(isapprox.(out.nodeflux[:],
    ComplexF64[15.190313925948457 + 0.0im, 3.0295193148601166e-6 - 1.8979297617655687e-8im, 6.835392097204295 + 0.0im, -2.3940362470180187e-6 + 1.4998131089887555e-8im, -6.835392097216894 + 0.0im, -5.116465430387518e-6 + 3.2053573997775784e-8im];
    rtol = 1e-8))

# output
true

JosephsonCircuits.hbnlsolve Method

julia

hbnlsolve(w::NTuple{N,Float64}, Nharmonics::NTuple{N,Int},
    sources::Vector{SourceTuple{N}}, psc::CompiledCircuit,
    circuitdefs::Dict{Any,Any}; kwargs...)

The general method on a compiled circuit psc, with the inputs in their canonical forms (tonefrequencies, sourcetable, definitiontable): it builds the mode set and the matrices at its mode count, then solves. It takes every keyword of the general method.

JosephsonCircuits.hbnlsolve Method

julia

hbnlsolve(w::NTuple{N,Number}, sources, frequencies::Frequencies{N},
    indices::FourierIndices{N}, psc::CompiledCircuit,
    nm::CircuitMatrices; kwargs...)

The nonlinear harmonic balance solve on an already compiled circuit psc with its matrices nm, at the mode set frequencies with its Fourier indices indices. This is what the other methods call after building those, with the inputs in their canonical forms (tonefrequencies, sourcetable); it takes every keyword of the general method except the ones which describe the mode set (Nharmonics, Nevaluationharmonics, maxintermodorder, frequencywindow, dc, odd, even) and the deprecated ones (ftol, switchofflinesearchtol, alphamin, factorization, maxharmonics and symfreqvar), and it does not accept method = Staged(), whose continuation builds each stage's own system. It takes one keyword the general method does not: reuse = nothing, an HBReuse which a solve fills and every later solve of the same topology, mode grid, tone frequencies, precision and backend rebinds (a mismatch is an ArgumentError); an ExternalSolver ignores it. It is refused with returnsystem = true and debugJacobian = true, whose system would be rebound under the caller.

Examples

julia

circuit = Circuit(
    [:p1 => Port(1; Z0 = :Rleft),
     :l1 => Inductor(:Lm),
     :l2 => Inductor(:Lm),
     :k1 => MutualInductor(:K1, :l1, :l2),
     :cc => Capacitor(:Cc),
     :jj3 => JosephsonJunction(:Lj),
     :jj4 => JosephsonJunction(:Lj),
     :cj => Capacitor(:Cj),
     :gnd => Ground()],
    [[(:p1, 1), (:l1, 1), (:cc, 1)],
     [(:cc, 2), (:l2, 1), (:jj4, 1), (:cj, 1)],
     [(:l2, 2), (:jj3, 1)],
     [(:p1, 2), (:l1, 2), (:jj3, 2), (:jj4, 2), (:cj, 2), (:gnd, 1)]])
circuitdefs = Dict{Symbol,Complex{Float64}}(
    :Lj =>2000e-12,
    :Lm =>10e-12,
    :Cc => 200.0e-15,
    :Cj => 900e-15,
    :Rleft => 50.0,
    :Rright => 50.0,
    :K1 => 0.9,
)

Idc = 50e-5
Ip=0.0001e-6
wp=2*pi*5e9
Nharmonics = (2,)
frequencies = JosephsonCircuits.removeconjfreqs(
    JosephsonCircuits.truncfreqs(
        JosephsonCircuits.calcfreqsrdft(Nharmonics),
        dc=true, odd=true, even=false, maxintermodorder=Inf,
    )
)
fi = JosephsonCircuits.fourierindices(frequencies)
Nmodes = length(frequencies.modes)
psc = JosephsonCircuits.compile(circuit)
nm = JosephsonCircuits.numericmatrices(psc, circuitdefs, Nmodes = Nmodes)

out=hbnlsolve(
    (wp,),
    [
        (mode=(0,),port=1,current=Idc),
        (mode=(1,),port=1,current=Ip),
    ],
    frequencies, fi, psc, nm)
# each entry, the direct current fluxes and the much smaller pump ones
all(isapprox.(out.nodeflux[:],
    ComplexF64[15.190313925948457 + 0.0im, 3.029519314860909e-6 - 1.8979297619056064e-8im, 6.835392097204296 + 0.0im, -2.394036246881364e-6 + 1.4998130573195207e-8im, -6.835392097216894 + 0.0im, -5.116465430251802e-6 + 3.2053573482853094e-8im];
    rtol = 1e-8))

# output
true

See the general hbnlsolve docstring for the formulation and the keywords.

JosephsonCircuits.hbsolve! Method

julia

hbsolve!(cache::HBCache, p::NamedTuple; warmstart = true)

Solve the nonlinear harmonic balance problem of cache at the design parameters p, warm starting from the last converged operating point; the parameters move the component values, the source currents and the tone frequencies written as their names. Returns the NonlinearHB solution; cache.converged reports whether it converged. A name of p which the definitions of the cache do not hold is an ArgumentError.

The compiled circuit and the mode grid are reused, and so are the system of the previous solve and what its method solves with, rebound to the new component values, drive and tone frequencies (see HBReuse): under NewtonKrylov the preconditioner and the Krylov vectors, under Newton and QuasiNewton the assembled Jacobian and its factorization, whose fill reducing ordering and symbolic analysis are kept, except for a circuit which injects direct current, whose canonical Jacobian Newton plans, orders and factorizes at every point. Otherwise only the numeric matrices and the solve itself are recomputed. The matrices are refilled on the patterns of the compiled circuit, which do not depend on the values, and assembled anew when a value changes the element type of its group (a resistance or a capacitance turned complex).

A start the parameters jumped away from can lie outside Newton's basin, so a warm started solve which does not converge is retried once from a cold start, the solve warmstart = false makes, keeping what the cache reuses. The warm attempt's messages reach the caller only when it is the outcome, so a retried point warns as its cold solve does, and a point costs at most the two solves, each within the solver's own budget. A retried point returns the cold solve, with the warm attempt's record ahead of its own in solverinfo.stages; where the circuit has several operating points, it is the one a cold start reaches, which need not be on the branch the sweep followed. A retry which converges becomes the stored point, and a point which fails both ways leaves the stored point as it was, so the next one starts from the last solution rather than from a non-solution or from nothing. warmstart = false starts cold, is not retried, and keeps the stored point, unlike JosephsonCircuits.reset!.

JosephsonCircuits.hbsolve Method

julia

hbsolve(ws, wp::NTuple{N,Any}, sources::Vector,
    Nmodulationharmonics::NTuple{M,Int}, Npumpharmonics::NTuple{N,Int},
    circuit, circuitdefs; dc = false, threewavemixing = false,
    fourwavemixing = true, maxpumpintermodorder = Inf,
    maxmodulationintermodorder = Inf,
    Nevaluationharmonics = map(i -> 2i, Npumpharmonics),
    frequencywindow = (0, Inf),
    maxmodulationharmonics = Nmodulationharmonics,
    iterations = 1000, atol = 1e-8, rtol = 0.0, method = NewtonKrylov(),
    x0 = nothing, warnnotconverged = true,
    nbatches = Base.Threads.nthreads(),
    returnS = true, returnSnoise = false,
    returnQE = true, returnCM = true, returnnodeflux = false,
    returnvoltage = false, returnnodefluxadjoint = false,
    returnvoltageadjoint = false, keyedarrays = true,
    temperature = 0.0, returnCnoise = false, returnVout = false,
    returnnbar = true, sensitivitynames::AbstractVector = String[],
    sensitivityoperatingpoint = true, sensitivitymode = :auto,
    returnSsensitivity = false, factorization = nothing,
    backend = CPU())

Solve a circuit driven by one or more strong pumps, then linearize about that operating point and sweep the weak signal frequencies ws. This calls hbnlsolve for the pump and hblinsolve for the signals and returns both solutions in an HB. Any number of pumps, ports and modes is supported, for three and four wave mixing.

The system is solved in a modified nodal analysis formulation in the node flux basis. Mutually coupled inductor branches are given auxiliary branch current variables, which keeps the system matrix bounded as a coupling coefficient approaches one, and one gauge fixing equation per floating inductive or Josephson subnetwork and zero frequency mode makes circuits with no inductive path to ground solvable without workaround inductors (see src/circuit/mna.jl). The nonlinear system is nondimensionalized by the scale Z0/w0 (see calcsolverscale), so atol does not depend on the unit system. The returned node fluxes and voltages contain only the node coordinates, not the auxiliary variables. The linearized solve throws an ArgumentError when any signal plus pump mode frequency is numerically zero; estimate a direct current limit from a sequence of decreasing nonzero frequencies instead.

Arguments

ws: the signal angular frequency or frequencies in radians per second, such as 2*pi*5.0e9 or 2*pi*(4.5:0.001:5.0)*1e9; a real number or any iterable of them.
wp::NTuple{N,Any}: the pump angular frequencies in radians per second, (2*pi*5.0e9,) for a single pump or (2*pi*5.0e9, 2*pi*6.0e9) for two, as real numbers or names, as a component value may be, which circuitdefs defines. The pumps should be non-commensurate; for commensurate pumps give the lowest frequency here and add the others to sources with a mode index equal to the frequency ratio. The pump solution's S then drives several port modes at once, so it is not a scattering matrix (see NonlinearHB); the linearized solution's is.
sources: the sources, as named tuples (mode, port, current) in any iterable. mode is a tuple of harmonic indices, one per pump frequency, port the port number, and current the complex current amplitude in Amperes; see sourcetable for the checks. [(mode=(1,0), port=1, current=Ip1), (mode=(0,1), port=1, current=Ip2)] applies a source at 1*wp[1] + 0*wp[2] and one at 0*wp[1] + 1*wp[2], both at port 1. A source with mode = (0,) and dc = true is a direct current bias. A current may be written in terms of parameters, as a component value is, and circuitdefs then gives them their values.
Nmodulationharmonics::NTuple{M,Int}: how many harmonics of each pump to retain around the signal in the linearized solve, which sets the signal and idler modes.
Npumpharmonics::NTuple{N,Int}: how many harmonics of each pump to retain as unknowns in the nonlinear solve. Its length is the number of non-commensurate pumps. The nonlinearity is evaluated on the larger Nevaluationharmonics grid.
circuit: a typed Circuit or a CompiledCircuit. A Circuit is compiled with the default node ordering; for another one, pass compile(circuit; sorting = ...).
circuitdefs: a dictionary from the symbols or symbolic variables used as component values or source currents to their numerical values. Optional when every component value and current is numeric.

Keywords

dc = false: retain the zero frequency mode in the nonlinear solve. A CurrentSource component of the netlist drives this mode with a constant current through itself from its first terminal to its second: it draws the current from the node at its first terminal and delivers it to the node at its second, the opposite sense of a port source, which injects its current into the port's first (positive) terminal. A nonzero one without the mode is an error.
threewavemixing = false: retain the even pump harmonics in the nonlinear solve, which are what three wave mixing processes couple through, and the signal modes offset from the signal by odd pump harmonics, the idlers of three wave mixing.
fourwavemixing = true: retain the odd pump harmonics in the nonlinear solve, and the signal modes offset by even pump harmonics, the idlers of four wave mixing.
maxpumpintermodorder = Inf: keep a pump mode which mixes two or more pumps only when its absolute harmonic indices sum to at most this order. Every harmonic of a single pump is kept up to its count in Npumpharmonics whatever the order, so with one pump it removes nothing (see truncfreqs).
maxmodulationintermodorder = Inf: the same truncation for the signal modes, every offset by the harmonics of a single pump kept up to its cap in maxmodulationharmonics.
Nevaluationharmonics = map(i -> 2i, Npumpharmonics): the harmonics of each pump on the grid where the nonlinearity is sampled, at least Npumpharmonics; twice the retained set by default, which dealiases the products of order three, the leading ones of a junction. See hbnlsolve.
maxpumpharmonics: deprecated and ignored with a warning; Npumpharmonics is the retained set and Nevaluationharmonics the sampling grid.
maxmodulationharmonics = Nmodulationharmonics: an upper bound on the absolute harmonic index retained for each pump in the signal modes, applied together with the intermodulation truncation; see truncfreqs.
frequencywindow = (0, Inf): a lower and upper bound, in the units of wp, on the absolute frequency of the retained pump modes; see hbnlsolve. The signal modes are not windowed.
iterations = 1000: the maximum number of nonlinear solver iterations before it returns unconverged.
atol = 1e-8: the residual tolerance norm(F) <= atol at which the nonlinear solution is considered converged. F is scaled by Z0/w0; see calcsolverscale. A tolerance below the rounding error of the scaled source, which a circuit whose interior is far from its port impedances can reach, is raised to it: no iteration resolves a residual below the rounding of the terms which make it.
rtol = 0.0: a relative residual tolerance; the nonlinear solve is converged when norm(F) <= max(atol, rtol*norm(F0)) with F0 the initial residual. It applies to the direct and Krylov methods; a Staged method ignores it, holding its final solve to atol and its interior stages to its own interioratol.
method = NewtonKrylov(): the nonlinear solver, an AbstractHBNonlinearSolver carrying its own options: NewtonKrylov (the matrix-free real Jacobian with a preconditioner, the default), Newton (the assembled real Jacobian), QuasiNewton (the complex holomorphic Jacobian, an approximation, with Anderson acceleration), Staged (source continuation on an adaptively grown harmonic grid, see stagedhbnlsolve, the strategy for operating points the direct methods fail outright; a search which ends without the point reports where it stalled, which does not prove that no operating point exists) or ExternalSolver.
x0 = nothing: an initial value for the node fluxes of the nonlinear solve, used by the direct and Krylov methods; a Staged method warm starts each stage from the last and refuses it.
warnnotconverged = true: warn when the nonlinear solve does not converge, and run the checks of a converged point which warn (a junction carrying nearly its critical current at direct current). A continuation whose stage solves are expected to fail passes false, and reports and checks its own outcome.
nbatches = Base.Threads.nthreads(): the number of batches the signal frequencies are split into for multithreading; 1 runs single threaded.
returnS = true: return the scattering parameters of the linearized solve.
returnSnoise = false: return the noise scattering parameters.
returnCnoise = false: return the added noise covariance at the output ports; see LinearizedHB.
returnVout = false: return the covariance of the waves leaving the ports, the inputs' noise at the temperatures of their terminations included.
returnnbar = true: return the occupation of each output port mode, the photons of the wave leaving it.
returnQE = true: return the quantum efficiency.
returnCM = true: return the commutation relations.
returnnodeflux = false, returnvoltage = false: return the node fluxes and voltages of the linearized solve.
returnnodefluxadjoint = false, returnvoltageadjoint = false: return the node fluxes and voltages of the adjoint linearized solve, the transposed system driven at each port and mode.
keyedarrays = true: return the outputs as keyed arrays with named, labeled axes rather than plain arrays.
temperature = 0.0: the physical temperature in kelvin of every dissipative element which does not state its own, and so of the noise it adds. A channel at temperature T carries the symmetrized noise nbar + 1/2, nbar its thermaloccupation, which at zero temperature is the vacuum's half photon. Raising it lowers the quantum efficiency and changes Cnoise but leaves Snoise and the commutation relations alone, since those describe the transformation rather than the state. It does not warm the ports: a port's termination states its own temperature, Port(n; termination = MatchedTermination(temperature = T)), zero unless given. A component may state its own temperature (Resistor(R; temperature = T), or a ScatteringParameters with noise = ThermalEquilibrium(T)), and takes this default otherwise. A block which states its noise with a NoiseCovariance, as an amplifier given by its scattering parameters does, adds that covariance whatever the temperature.
sensitivitynames = String[]: the components to take sensitivities with respect to, named by their identifiers as symbols or strings. Supported types are C, L, R and Lj with numeric values. A ScatteringParameters block has no scalar value to perturb and cannot be named; see designsensitivities for sensitivities with respect to the parameters of a block.
sensitivitypairs, sensitivityblockpairs, nsensitivityparameters, sensitivitylabels: the design parameter interface used by designsensitivities, which names physical parameters rather than components; each pair (componentname, parameterindex, alpha) gives the relative direction alpha = (dv/dp)/v of a component value under a parameter. Not intended to be passed directly.
sensitivitymode = :auto: the order in which the operating point shift is contracted into the sensitivities. :forward costs one product against the linearized system per component and signal frequency; :reverse pushes the output functionals through the transposed pump Jacobian once per output port and mode pair, so its cost does not grow with the number of components. :auto chooses :reverse when the components outnumber the output port and mode pairs eight times over, or when the forward order's stamps, a value per stored entry of the linearized system per component, would exceed FORWARDSENSITIVITYSTAMPBYTES. Both take the parameters of scattering blocks as they take components, and support any number of pumps.
sensitivityoperatingpoint = true: include the shift of the pump operating point in the sensitivities, making them total derivatives rather than derivatives at a fixed operating point. Near the gain peak of a strongly pumped amplifier the operating point term is comparable to or larger than the fixed point term. Requires the exact real Jacobian of the nonlinear solution, which is assembled once and factorized once for the forward contraction order and once per batch of signal frequencies for the reverse one, whose transposed solves run on the batch's thread: a sparse factorization cannot be solved against from several threads at once. Without Josephson junctions the operating point contribution is identically zero and is skipped.
returnSsensitivity = false: return dS/dr, the derivative of the scattering matrix with respect to a relative perturbation r of each named component value (p -> r*p at r = 1), computed by the adjoint method.
factorization = nothing: the factorization of the linearized solve at each signal frequency. nothing takes the backend's sparse factorization for one tone, KLUfactorization on the host and cuDSS (CUDSSFactorization) on a device, and for two or more tones BlockFactorization, the dense node blocks of the multi-tone system with BLAS-3 arithmetic, when its factors fit in half the free memory (linearizedfactorization). Any of them can be given explicitly. The block factorization pivots only within its dense blocks and stops at a singular one: chosen by nothing, the sweep then runs again on the sparse factorization, and given explicitly it throws. On a device the choice also picks the solver of the batch, cuDSS for a sparse factorization and the batched block factorization for a BlockFactorization. QRfactorization is refused (see hblinsolve). The nonlinear solve's factorization is an option of its method. The precision of the linearized solutions is the factorization's: BlockFactorization(precision = Float32, refine = false) solves each signal frequency entirely in single precision (see BlockFactorization). The nonlinear solve's precision is NewtonKrylov(precision = ...).
backend = CPU(): the KernelAbstractions backend both solves run on. The nonlinear solve assembles, factorizes and iterates there; the linearized sweep solves batches of signal frequencies there and falls back to the host for what it cannot serve (see hblinsolve).
ftol: deprecated, read as atol with a warning.
switchofflinesearchtol, alphamin: deprecated and ignored with a warning.
symfreqvar = nothing: deprecated, the parameter a frequency dependent value was written as an expression in. Write the value as a FrequencyDependent closure of the frequency instead.
returnZ, returnZadjoint, returnZsensitivity, returnZsensitivityadjoint: removed; passing any of them warns. Compute impedances from the scattering parameters instead.

Returns

HB: the nonlinear and linearized solutions; see HB.

JosephsonCircuits.hbsolve Method

julia

hbsolve(ws::Vector{Float64}, wp::NTuple{N,Float64},
    sources::Vector{SourceTuple{N}}, Nmodulationharmonics::NTuple{M,Int},
    Npumpharmonics::NTuple{N,Int}, psc::CompiledCircuit,
    circuitdefs::Dict{Any,Any}; kwargs...)

The method on a compiled circuit psc with the inputs in their canonical forms (sweepfrequencies, tonefrequencies, sourcetable, definitiontable). The method on a Circuit converts its inputs and calls this one, so that the solve is compiled once for every way of writing them; it takes the keywords documented for that method.

JosephsonCircuits.hbstability Function

julia

hbstability(circuit, circuitdefs = Dict(); nonlinear = nothing,
    Nmodulationharmonics = (8,), method = Monodromy(),
    pumpfrequency = nothing, frequencyscale = nothing,
    factorization = KLUfactorization())

The local stability of the periodic operating point nonlinear, a converged hbnlsolve solution with one pump frequency, or of zero junction phase without one: the temporal poles of the circuit linearized about it, as an HBStabilityResult, whose real parts tell whether a small perturbation of the orbit grows. They say nothing of what larger perturbations do, and settle the orbit's stability only as far as the truncation, in harmonics or steps, is converged. The circuit and its definitions must be those of the operating point; commensurate drives are harmonics of its one pump.

The perturbation is expanded in every harmonic -H:H of Nmodulationharmonics = (H,), whatever the pump's symmetry, over the junction modulation hblinsolve samples from the same operating point; without an operating point or pumpfrequency it has the one harmonic 0. pumpfrequency gives the fundamental of pumped blocks without an operating point, and must agree with one. The polynomial methods truncate the perturbation at these harmonics: converge their poles in the pump's harmonics and in these, which edgeweights and a solve at a larger H tell. The period map, Monodromy, takes them as the width of its profiles alone, requiring 2H + 1 <= steps. Converge its poles in the pump's harmonics and in its steps. Its rateerrors is an estimate from the map at twice the steps, not a bound; near zero growth, compare full solves at increasing steps.

Choose a method for the desired spectral coverage:

Monodromy() returns up to 10 least damped resolved modes of a pumped circuit, one representative each, including line delays. Set nev to change that count. Without a pump it uses DenseSpectrum() with the default 1000-unknown limit.
DenseSpectrum returns accepted finite poles of a small harmonic truncation, including aliases.
ShiftInvert searches near specified complex frequencies.
ContourIntegral searches inside a circle. It is required for LaplaceResponse models and exact delays without a pump.

A shift or contour search can miss instabilities outside its region. factorization factors the polynomial on the host, KLUfactorization or LUfactorization, refactoring its fixed pattern, and the transient's steps for the period map on the host. frequencyscale, an angular frequency, scales the polynomial: the pump frequency by default, or a natural frequency of the unpumped circuit; it changes no pole.

The variational equation is Q(s) = K + G (s + D) + C (s + D)^2, with D = diag(im*n*pumpfrequency) and the junction stiffness along the orbit in K. Floating flux references are replaced by voltages to retain physical zero poles and RC decay. Ports keep their terminations, and independent sources are held fixed.

Supported models include real constant lumped elements, mutual inductors, junctions, nonlinear inductors, constant and RationalScattering blocks with internal states, and LinearizedScattering blocks. Use a polynomial method for constant/rational conversion harmonics, or ContourIntegral for analytic Laplace harmonics. Pumped Monodromy requires a supported time-domain realization. Tabulated data or a callable of real frequency needs a causal rational fit or LaplaceResponse; complex constant lumped values are refused. Refine a positive pole in harmonics or steps before interpreting it as an instability of the physical circuit.

Examples

julia

julia> c = Circuit([(:r, 1, 0, Resistor(2.0)),
                   (:c, 1, 0, Capacitor(0.5)),
                   (:l, 1, 0, Inductor(2.0))]);

julia> p = hbstability(c);

julia> sort(imag.(p.poles)) ≈ [-sqrt(3)/2, sqrt(3)/2]
true

julia> all(isapprox.(real.(p.poles), -0.5)) && p.converged
true

JosephsonCircuits.matchpoles Method

julia

JosephsonCircuits.matchpoles(reference, candidates;
    indices = eachindex(reference.poles),
    candidateindices = eachindex(candidates.poles),
    minoverlap = 0.9, maxdistance = Inf, aliases = true)

Match selected reference poles to distinct candidates by maximizing the total physical profile overlap, using node voltages and junction fluxes. Works across different harmonic truncations: absent harmonics are zero padded, with their tails retained in the normalization. Node and junction axes must agree. Voltages on both sides use the reference frequency scale.

Returns one named tuple per selected reference pole: index (zero when unmatched), overlap, alias and distance. With aliases = true and equal nonzero pump frequencies, the candidate's nearest Floquet alias is compared: candidate.poles[index] - im*alias*pumpfrequency. Its harmonic coefficients are shifted by alias at the same time. distance is the absolute difference of these physical rates, in inverse seconds.

minoverlap rejects weak matches and maxdistance optionally limits rate motion. candidateindices selects acceptable candidates without copying a result; returned indices still refer to the original result. Unmatched poles are reported explicitly; a parameter continuation should reduce its step or expand its search. An overlap is not a proof of branch identity at a degeneracy. Remove duplicate discoveries and boundary artifacts from the candidate set before interpreting matches. Modes with no observable node voltage or junction flux cannot be matched this way. No new circuit solve is launched by this function.

Examples

julia

julia> c = Circuit([(:r, 1, 0, Resistor(2.0)), (:c, 1, 0, Capacitor(0.5))]);

julia> p = hbstability(c);

julia> matched = only(JosephsonCircuits.matchpoles(p, p; indices = [1]));

julia> matched.index == 1 && matched.overlap ≈ 1 && matched.distance == 0
true

JosephsonCircuits.reset! Method

julia

reset!(cache::HBCache)

Discard the stored operating point, so the next hbsolve! starts cold, and with it what the previous solves taught the preconditioner: a coupling set grown by escalation or by measurement, and the deflation candidates of a Floquet preconditioner, so the next solve builds the preconditioner its method asks for. A warm start which fails is retried cold without discarding any of this (see hbsolve!); use it when the parameters jump far enough that the previous solution is a worse starting point than zero, which saves that failed attempt, when crossing to a different solution branch, or to drop what the preconditioner grew.