Fixed schedules and repeatable propagation#

A fixed schedule prescribes where the solver steps. Use it to replay positions or run compatible problems in a batch, after checking its accuracy for the parameter range. Accepting every prescribed step does not make its local error acceptable. Adaptive stepping remains the starting point for discovering a suitable schedule.

Fixed schedules: discover, calibrate, replay#

Adaptive stepping does not remove the convergence study, and fixed stepping is not less accurate. It wins outright in three places: for ensembles, where identical steps make the numerical error a common systematic that cancels between members instead of a per-member variance that adds; for repeat calculations and debugging, since a replayed run takes the recorded steps and meets the original run’s accuracy target (it need not reproduce its bits, SEP 25); and on a device, where there is nothing to synchronise. Both modes are therefore first class, and batch lockstep execution uses fixed schedules only.

Discover. One adaptive run records \(h(z)\); discover turns the accepted step lengths into a FixedSchedule.
Calibrate. calibrate replays it at \(h\), \(h/2\) and \(h/4\) and reads the observed order from the ratio of successive differences, then estimates the error the base schedule actually carried by Richardson extrapolation.
Replay. replay rescales the schedule along the measured exponent to reach a target error. Rescaling follows the discovered profile, so wherever the adaptive run chose to go slowly, the replayed one does too.

>>> import numpy as np
>>> from solprop.stepper import FixedSchedule, calibrate, replay
>>> base = FixedSchedule(np.full(8, 0.125))
>>> # Stand-in for a run: an observable whose error falls as h**5.
>>> quantity = lambda schedule: 1.0 + schedule.max_step**5
>>> calibration = calibrate(base, quantity, order=5, requested=1e-8)
>>> round(calibration.observed_order, 6)
5.0
>>> replay(base, calibration, calibration.achieved_error / 32).n_steps
16

Two refusals are built into this. The order gate rejects an observed order more than three above the scheme’s: refining removed something that is not truncation error, so the coarse schedule is off its convergence curve rather than merely large. Three is the margin because a legitimate strongly nonlinear stretch reads about half an order high while this failure mode reads nine orders high, and a two-level extrapolant cannot tell them apart at all. Separately, when the coarsest difference sits within a factor of 1000 of the round-off floor, calibrate marks the result inconclusive and declines to pronounce an order — below the floor the differences measure accumulated rounding, and the same ratio has come back as \(-0.83\) there and been read as healthy. An inconclusive calibration cannot be replayed from; that is the point of measuring rather than assuming.

Replaying a run’s own steps#

Two functions turn an adaptive run into a schedule. discover takes the accepted step lengths as a list. FixedSchedule.from_history reads the step_history a run writes to its output, which records every attempt with its start, its length and whether it was accepted, and drops the rejected attempts. Replaying the schedule with the same save planes repeats the run’s step positions without its controller. With a method whose orders are fixed (below), that is how two precisions or two backends are compared on one discretisation; prop_gnlse steps with RungeKutta() (a few seconds):

import numpy as np

from solprop.driver import prop_gnlse
from solprop.stepper import FixedSchedule


def soliton(**changes):
    """Run a third-order soliton over 0.2 m of fibre."""
    return prop_gnlse(
        0.11,
        0.2,
        [0.0, 0.0, -11.83e-27],
        wavelength=835e-9,
        wavelength_limits=(500e-9, 2e-6),
        trange=2e-12,
        width=100e-15,
        peak_power=9 * 11.83e-27 / (0.11 * 1e-26),
        shape="sech",
        n_save=3,
        noise=None,
        progress="none",
        # In both runs: the default in-stage collar would refuse this history (below).
        absorber_scheme="split_step",
        **changes,
    )


def steps_of(run):
    """Return a run's accepted steps as a schedule."""
    return FixedSchedule.from_history(
        np.asarray(run["step_history/z"]),
        np.asarray(run["step_history/h"]),
        np.asarray(run["step_history/accepted"]),
    )


schedule = steps_of(soliton())
one = soliton(schedule=schedule)
other = soliton(schedule=schedule)  # change backend= or precision= here to compare
print(np.array_equal(steps_of(one).steps, schedule.steps))  # True

What a replay repeats. Under RungeKutta() a step is fixed by its length, so a replay takes the adaptive run’s discretisation, and so does LowStorageRungeKutta(). The default Multistep() chooses its order as it goes: after each accepted step it weighs the neighbouring orders’ defect estimates and adopts one only when it beats the current order by a hysteresis margin. A replay keeps the schedule’s step positions but makes that choice afresh, so its orders can differ in two ways:

from the adaptive run’s, where that run rejected a step: it then waits a step before raising the order, which a replay never has to;
from another replay’s in another precision or on another backend, where an estimate lies near the margin: rounding decides the choice there, and a different order sequence makes a different truncation error.

Two replays in one precision on one backend are identical. To compare the arithmetic alone, replay with fixed orders: RungeKutta(), as above, or Multistep(variable_order=False), whose order rises from one as its history accrues to order and stays there, by a rule that reads no estimate. Two such replays of one schedule make the same truncation errors, so their difference is rounding alone. Single precision steps with RungeKutta(), and refuses Multistep(), so a comparison across precisions is a Runge–Kutta one. Recording the order sequence in the schedule would also fix a variable-order replay’s orders, and make it match its adaptive run; that would extend the saved schedule and is not done, since a fixed order gives the comparison without it.

The collars. Under the default absorber_scheme="nonlinear" the temporal collar’s largest rate is an eigenvalue of the step, so a replayed step may be no longer than the collar’s stable step (stable_step, scaled to the method’s stability interval), and translation refuses a schedule that has one, naming the limit. An adaptive run is not held to that limit, because an empty collar produces no error, so a default run’s own history is often refused: on the soliton above its longest accepted step is 5.60 mm against a limit of 5.24 mm under RungeKutta(), and 1.59 mm against 0.139 mm under Multistep(). Either of two remedies replays exactly:

give absorber_scheme="split_step" to both the discovery run and the replay, as above (on the replay alone the replay runs, but differs from the discovery run by the change of scheme); or
hold the discovery run to a dz_max at or under the limit the refusal prints, and replay its history with the same dz_max. The printed limit is rounded down, so the printed value itself is safe.

A stability limit belongs to the integrator, Multistep()’s being far shorter than RungeKutta()’s, so run the discovery with the integrator the replay will use: the printed limit is computed for it.

A requested spatial collar on a spatial route (SEP 40) adds three more limits on a replayed step: its own stable step, which also depends on the integrator, and the fixed-schedule crossing and damping bounds, which do not. An adaptive run is held only to the looser adaptive crossing limit, so its history can break them. The spatial collar is always in-stage, so absorber_scheme="split_step" does not lift these; dz_max at or under the printed limit in the discovery run does, and each refusal says so.

Two more limits. A replay lands on its save planes, so interpolate_planes=True is refused beside a schedule; and a discovery run with interpolate_planes=True steps past its planes, so its history does not contain them and is refused as a schedule for the same planes. The accepted steps must also follow each other, each starting where the last ended. A resumed run’s history does when the resumed part continues from where the history ends, as it does after a stop at a save plane; a history that does not, such as two runs’ joined, is refused.

Plain fixed steps#

A fixed-step run with no adaptivity at all needs no discovery. uniform cuts the run at every save plane and hard kink and covers each segment with the fewest equal steps of at most dz, so no step exceeds dz and every target is an edge of the schedule. There is no dense output, so a fixed schedule must contain the planes; this constructor guarantees it.

>>> from solprop.stepper import FixedSchedule
>>> FixedSchedule.uniform(0.0, 1.0, 0.3, save_planes=[0.5]).steps
array([0.25, 0.25, 0.25, 0.25])

Replaying a schedule is not blind. The embedded error estimate costs nothing extra, so a fixed run still computes the scaled norm every step, records it in the step history and summarises it in ControlRecord; when met_tolerance is false a step was too coarse for the tolerance the norm is scaled by, and the driver warns. That is the fixed-mode counterpart of the dz_min error: neither mode fails silently.

Ensembles with noise. Calibrate a schedule for a noisy ensemble with noise on, and check it on more than the seed that discovered it. A schedule discovered on a noiseless run integrates the seed’s own dynamics poorly: on a 20 pJ fibre soliton it got the noise-driven part of the field wrong by 6 % with RungeKutta() and 78 % with Multistep() (SEP 37, §3.4).

The low-storage march#

integrator=LowStorageRungeKutta() (choosing a stepper) replays a schedule and does nothing else: it has no embedded estimate, so there is no norm to report and none to warn on. Every step’s recorded norm is 0.0, meaning “not estimated”, and met_tolerance is then vacuously true. Its accuracy is calibrated from outside, on the quantity you will report, with calibrate(..., order=4):

Calibrate the observable, not rtol. The march’s rtol sets nothing but the recorded tolerance and the spatial collar’s diagnostic threshold. The observable can be anything the run produces: an energy in a band orders below the pump, a peak, a difference of two strong fields. A whole-field norm sees a band holding a fraction \(\delta\) of the energy at about \(\sqrt\delta\) of its weight, which is why calibrating on the band itself is the safer check.
requested is a label. A uniform schedule was not discovered by an adaptive run, so there is no tolerance it was asked for; pass the controller’s rtol or the error you are aiming at, which calibrate only quotes in its message.
The order gate reads four. With order=4 a ladder observing more than seven is refused: refining removed something that is not truncation error. In the asymptotic range (\(|\operatorname{Im} D|\,h/2 \lesssim 1\) in the populated bins) the march reads four; a coarser base reads lower and is not yet a calibration.
A weak observable needs its floor. The strong field’s rounding leaks into a weak band and sets the floor calibrate(floor=...) must be told. The difference between two backends on one schedule is pure rounding and measures it: native against numpy_ref, or one backend’s two flow treatments (propagator_mode="tabulate" and "recompute").
The flow table follows the schedule. Tabulating a constant operator holds one copy of the state, \(e^{\bar h D/2}\), built for the mean step \(\bar h\) of each run of equal steps and reused along it: a uniform segment builds it once, and the phase it applies over the segment is the segment’s own. A schedule whose steps all differ (a discovered one) rebuilds it every step, where propagator_mode="recompute" is the cheaper choice.

A third-order soliton’s first compression, scored on the light below 780 nm, which holds less than a millionth of the energy (a few seconds):

import numpy as np

from solprop.driver import Simulation
from solprop.stepper import FixedSchedule, LowStorageRungeKutta, calibrate, replay

length = 0.2  # m


def weak_band(schedule):
    """Replay `schedule` and return the energy fraction below 780 nm."""
    simulation = Simulation.from_gnlse(
        0.11,
        length,
        [0.0, 0.0, -11.83e-27],
        wavelength=835e-9,
        wavelength_limits=(500e-9, 2e-6),
        trange=2e-12,
        width=100e-15,
        peak_power=9 * 11.83e-27 / (0.11 * 1e-26),
        shape="sech",
        n_save=2,
        noise=None,  # a noiseless soliton: the calibration is about the stepper
        schedule=schedule,
        integrator=LowStorageRungeKutta(),
        # The default in-stage collar caps an explicit step; the split one does not.
        absorber_scheme="split_step",
        dz_max=schedule.max_step,
        dz_init=schedule.max_step,
    )
    output = simulation.run(progress="none")
    spectrum = np.abs(output["field"][-1, 0, 0]) ** 2
    wavelength = 2 * np.pi * 299792458.0 / simulation.grid.omega
    return float(spectrum[wavelength < 780e-9].sum() / spectrum.sum())


base = FixedSchedule.uniform(0.0, length, length / 32)
calibration = calibrate(base, weak_band, order=4, requested=1e-6)
print(f"observed order {calibration.observed_order:.2f}")  # 3.92
production = replay(base, calibration, calibration.achieved_error / 10)
print(production.n_steps)  # 58

A replay is accepted in place, so an exception or a signal inside a step leaves no state to checkpoint: the run says so and its file ends at the last saved plane, from which resume=True continues; when that plane kept no field, the state held there is checkpointed instead (field retention). Save planes often enough (step_planes=k saves the end of every k-th step) to bound what a stop loses.