solprop.backends.exponential#

The phi functions and exponential Adams weights, in NumPy (SEP 16 sections 4.2-4.3).

Exponential Adams integrates the polynomial through untransported source samples exactly against the linear flow. For a constant operator D, a step h and the Lagrange basis ell_j(theta) = sum_k c_jk theta**k on the normalised nodes,

w_j(x) = int_0^1 exp((1 - theta) x) ell_j(theta) dtheta
= sum_k c_jk k! phi_{k+1}(x),

with x = h D per element and phi_r(x) = sum_m x**m / (m + r)!. At a fraction xi of the step the same integral to xi is sum_k c_jk k! xi**(k+1) phi_{k+1}(xi x). At x = 0 these are the ordinary Adams weights (solprop.stepper.multistep.integrated_weights()).

The phi functions are evaluated by two algorithms, each where it is stable:

near zero (|x| <= 8): scale to y = x / 2**s with |y| <= 1/2, sum the Taylor series of the highest function to a fixed number of terms, recur down with phi_k(y) = 1/k! + y phi_{k+1}(y), which contracts an error by |y| and never divides by y, and undo the scaling by the doubling identity phi_k(2y) = [exp(y) phi_k(y) + sum_{j=1}^k phi_j(y) / (k - j)!] / 2**k, obtained by splitting the defining integral at one half [1] [2];
away from zero (|x| > 8): recur up from phi_0 = exp(x), phi_k(x) = (phi_{k-1}(x) - 1/(k-1)!) / x, whose division by |x| > 8 contracts the error of every earlier function.

Both are checked against augmented matrix exponentials and independent quadrature in tests/backends/test_exponential_weights.py, on both sides of the join. This module is the reference the compiled backends’ moment kernels are checked against.

References#

[1]

S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems”, J. Comput. Phys. 176, 430 (2002). doi:10.1006/jcph.2002.6995

[2]

M. Hochbruck and A. Ostermann, “Exponential integrators”, Acta Numerica 19, 209 (2010), section 2. doi:10.1017/S0962492910000048

Module Attributes

`MAX_PHI_ORDER`	Highest phi function evaluated: the corrector's `MAX_ADAMS_ORDER + 1` nodes need `phi_1` to `phi_{MAX_ADAMS_ORDER + 1}`, and one more serves a neighbouring order.
`PHI_JOIN`	`\|x\|` above which the phi functions are recurred up from `exp(x)`.
`COUPLED_SECTIONS`	The coupled flow's section methods the exponential formulation steps (SEP 16 Q3).
`EXPONENTIAL_REPRESENTATIONS`	The linear-operator representations the exponential formulation steps (SEP 16).

Functions

`carrier_chord`(linear, index)	Return the reference element's operator alone, or `None` for a separable one.
`carrier_reference`(problem, method)	Return where a run's carrier frame takes its shift, or `None` when off.
`exponential_weights`(nodes, x[, upper, helpers])	Return the exponential Adams weights of distinct nodes at each argument.
`phi_functions`(x, order)	Return `phi_0(x)` to `phi_order(x)` for each element of x.
`refuse_carrier_frame`(problem, backend, *[, ...])	Refuse a carrier frame where it is not yet delivered (SEP 35 section 5).
`refuse_formulation`(problem, backend)	Refuse the exponential formulation where it is not yet delivered (SEP 16).

Classes

`CarrierChord`([constant, expansion, table])	The reference element's cumulative operator alone (SEP 35 sections 6.1, 6.4).
`CarrierReference`(policy, bin, omega, reason)	Where a run's carrier frame takes its shift, fixed for the run (SEP 35 6.1).

solprop.backends.exponential.COUPLED_SECTIONS = frozenset({'constant_schur', 'frozen_residual'})#: The coupled flow’s section methods the exponential formulation steps (SEP 16 Q3).

class solprop.backends.exponential.CarrierChord(constant=None, expansion=None, table=None)[source]#

Bases: object

The reference element’s cumulative operator alone (SEP 35 sections 6.1, 6.4).

The operator’s own representation restricted to column 0, polarisation 0 and the reference bin, every batch member kept: its chord over a step is the reference element’s exponent, to rounding the one the backend forms, at one scalar per member. The step cap reads it for every proposal, and a backend that forms its exponent elsewhere (the native engine, on a device too) takes its shift from it with nothing read back.

Attributes:

constantnumpy.ndarray or None: D at the element, 1/m, shape (n_batch_c,), for a constant operator.
expansionsolprop.translate.expansion.ParameterExpansion or None: The expansion with its terms at the element, for an expansion.
tablesolprop.translate.chebyshev.PiecewiseChebyshev or None: The z table at the element, for a table.

constant: Complex[ndarray, 'n_batch_c'] | None = None#

expansion: ParameterExpansion | None = None#

exponent(z1, z2)[source]#

Return Phi(z2) - Phi(z1) at the element, per batch member.

Parameters:

z1, z2float: Positions in m, inside the operator’s span.

Returns:

numpy.ndarray: Dimensionless, shape (n_batch_c,).

table: PiecewiseChebyshev | None = None#

class solprop.backends.exponential.CarrierReference(policy, bin, omega, reason, chord=None)[source]#

Bases: object

Where a run’s carrier frame takes its shift, fixed for the run (SEP 35 6.1).

The shift of batch member b is sigma_b = i Im x[b, 0, 0, bin], x the trial’s own exponent Phi(z + h) - Phi(z): column 0 is the frame’s mode (build_phi fixes the frame on the first mode; a coupled route’s carrier is the same in every channel), polarisation 0. The change of variable is exact for any shift, so the choice of bin only changes how smooth the interpolated source is.

Attributes:

policy{‘scalar’}: The policy asked for (solprop.stepper.multistep.CarrierFrame).
binint or None: The active bin nearest the reference frequency, ties to the lower index; None where the policy is inert (sigma = 0 for the run).
omegafloat or None: That bin’s angular frequency, rad/s; None with bin.
reasonstr: How the reference was chosen, or why the policy is inert, for the run’s log.
chordCarrierChord or None, optional: The reference element’s operator alone (carrier_chord()), for the step cap and a backend that forms its shift on the host; None when inert.

property active: bool[source]#

Return whether the shift can be nonzero.

Returns:

bool: False when the run takes sigma = 0, the present arithmetic exactly.

bin: int | None#

chord: CarrierChord | None = None#

chord_phase(z1, z2)[source]#

Return the largest carrier phase over the batch from the chord alone.

Parameters:

z1, z2float: The proposed step’s ends, m.

Returns:

float: max_b |Im (Phi(z2) - Phi(z1))| at the reference element, rad; 0 when the policy is inert.

Raises:

BackendError: If the reference has no chord (a separable operator).

omega: float | None#

phase(exponent)[source]#

Return the largest carrier phase over the batch, for the step cap.

Parameters:

exponentnumpy.ndarray: x over a proposed step.

Returns:

float: max_b |Im x[b, 0, 0, bin]|, rad; 0 when the policy is inert.

policy: Literal['off', 'scalar']#

reason: str#

shift(exponent)[source]#

Return each batch member’s shift for a trial’s exponent.

Parameters:

exponentnumpy.ndarray: x = Phi(z + h) - Phi(z), in the canonical layout.

Returns:

numpy.ndarray: Purely imaginary, shape (n_batch, 1, 1, 1), broadcasting against the state.

Raises:

BackendError: If the policy is inert: such a run takes the unshifted arithmetic.

solprop.backends.exponential.EXPONENTIAL_REPRESENTATIONS = ('constant', 'separable', 'expansion', 'table')#

The linear-operator representations the exponential formulation steps (SEP 16).

A constant operator’s reference flow over a step is its exact flow and its residual is the source alone, R_n = F. A varying one (expansion, table) takes the compensated reference of SEP 16 section 4.2 (phase 4): the chord over the step and the transported compensation (D(z_j) - B_n) u_j, with D the derivative of the representation itself.

solprop.backends.exponential.MAX_PHI_ORDER = 8#: Highest phi function evaluated: the corrector’s MAX_ADAMS_ORDER + 1 nodes need phi_1 to phi_{MAX_ADAMS_ORDER + 1}, and one more serves a neighbouring order.

solprop.backends.exponential.PHI_JOIN = 3.0#

|x| above which the phi functions are recurred up from exp(x).

Below it the upward recurrence loses digits to the subtraction of 1/(k-1)!: at order eight, 291 ulp of the natural scale at |x| = 2 and 32 at 3, against augmented exponentials. Above it the scaled series would need more than three doublings. Three keeps both sides within the witnesses’ 256 ulp with a margin and bounds the doublings, the most expensive part of the near branch (SEP 16 phase 3; eight, the SEP 15 study’s join, cost up to four).

solprop.backends.exponential.carrier_chord(linear, index)[source]#

Return the reference element’s operator alone, or None for a separable one.

Parameters:

linearsolprop.translate.linop.LinearOperator: The cumulative operator.
indexint: The reference bin.

Returns:

CarrierChord or None: None for the separable (spatial) representation, which the carrier frame refuses.

solprop.backends.exponential.carrier_reference(problem, method)[source]#

Return where a run’s carrier frame takes its shift, or None when off.

The reference rule of SEP 35 section 6.1. An envelope grid without the third harmonic subtracts beta_0, so its shift would be the table’s rounding: the policy is inert there and the run is the present arithmetic exactly. An envelope grid with the third harmonic keeps the carrier phase whole, its operator is the field grid’s, and it is treated as one. The reference frequency is the grids’ omega0; a custom frame with another reference changes only how much of the rotation the shift removes, never the step’s exactness.

Parameters:

problemsolprop.translate.problem.KernelProblem: The problem.
methodsolprop.stepper.multistep.Multistep: The method.

Returns:

CarrierReference or None: None for carrier_frame="off" or the interaction formulation.

solprop.backends.exponential.exponential_weights(nodes, x, upper=1.0, *, helpers=Helpers(name='numba'))[source]#

Return the exponential Adams weights of distinct nodes at each argument.

w_j = int_0^upper exp((upper - theta) x) ell_j(theta) dtheta, so that exp(upper x) u_n + h sum_j w_j F_j advances u' = D u + F by upper h when F is the polynomial interpolating the samples F_j at the nodes and x = h D.

Parameters:

nodesnumpy.ndarray or list of float: Distinct normalised positions (z_j - z_n) / h.
xnumpy.ndarray or complex: h D per element, dimensionless.
upperfloat, optional: The fraction of the step integrated, in [0, 1].
helperssolprop._helpers.Helpers, optional: The implementation of the Lagrange basis (SEP 29 section 6.6); the shared Numba one by default.

Returns:

numpy.ndarray: One weight per node along a new leading axis, each of x’s shape.

Raises:

StepperError: If the nodes are not distinct and finite, upper lies outside [0, 1] or an argument overflows.

Examples

At x = 0 the weights are the ordinary Adams-Moulton ones:

>>> import numpy as np
>>> from solprop.backends.exponential import exponential_weights
>>> (24 * exponential_weights([1, 0, -1, -2], 0.0)).real.round(12).tolist()
[9.0, 19.0, -5.0, 1.0]

solprop.backends.exponential.phi_functions(x, order)[source]#

Return phi_0(x) to phi_order(x) for each element of x.

phi_0(x) = exp(x) and phi_r(x) = int_0^1 exp((1 - theta) x) theta**(r-1) / (r-1)! dtheta.

Parameters:

xnumpy.ndarray or complex: The arguments, dimensionless (a step times the linear operator). Their real parts must not exceed 700, where exp overflows.
orderint: The highest function, in [1, MAX_PHI_ORDER].

Returns:

numpy.ndarray: phi_r(x) along a new leading axis of length order + 1.

Raises:

StepperError: If an argument is not finite or overflows, or the order is out of range.

Examples

>>> import numpy as np
>>> from solprop.backends.exponential import phi_functions
>>> values = phi_functions(np.array([0.0, 1.0]), 2)
>>> np.round(values.real, 12).tolist()
[[1.0, 2.718281828459], [1.0, 1.718281828459], [0.5, 0.718281828459]]

solprop.backends.exponential.refuse_carrier_frame(problem, backend, *, delivered=frozenset({}))[source]#

Refuse a carrier frame where it is not yet delivered (SEP 35 section 5).

Parameters:

problemsolprop.translate.problem.KernelProblem: The problem.
backendstr: The backend’s name, for the message.
deliveredfrozenset of str, optional: The policies this backend steps; none by default.

Raises:

BackendError: Naming the policy and what does not deliver it.

solprop.backends.exponential.refuse_formulation(problem, backend)[source]#

Refuse the exponential formulation where it is not yet delivered (SEP 16).

The exponential formulation steps constant, separable and varying diagonal operators (phases 2–4) and the constant-Schur and frozen-residual sections of a coupled flow (phase 5); a CF4 section is excluded (Q3) and refused by name before any source is evaluated.

Parameters:

problemsolprop.translate.problem.KernelProblem: The problem.
backendstr: The backend’s name, for the message.

Raises:

BackendError: Naming what is unsupported.