diff --git a/changelog.d/per-band-loss-interval.fixed.md b/changelog.d/per-band-loss-interval.fixed.md new file mode 100644 index 000000000..41729ef21 --- /dev/null +++ b/changelog.d/per-band-loss-interval.fixed.md @@ -0,0 +1 @@ +The boundaries guide reads a band's absorbed energy from saved planes with an exact interval formula in place of the trapezoid, which overcounted a strongly absorbed band's loss by up to twice; the medium's own absorption is now its own share (#573). diff --git a/changelog.d/walkoff-periodic-sum.fixed.md b/changelog.d/walkoff-periodic-sum.fixed.md new file mode 100644 index 000000000..08196cc89 --- /dev/null +++ b/changelog.d/walkoff-periodic-sum.fixed.md @@ -0,0 +1 @@ +The temporal collar's walk-off integral, read by `warn_walkoff` and by `absorber(grid, walkoff=...)`, is now the rectangle sum over the periodic time axis, which is exact for the taper, rather than a trapezoid that left out one sample; a collar sized for a drift now removes exactly the stated `exp(-2Γ)` of its power, and the warning fires at the exact threshold (#490). diff --git a/docs/dev/kernel-contract.md b/docs/dev/kernel-contract.md index c9c464d05..40a0779be 100644 --- a/docs/dev/kernel-contract.md +++ b/docs/dev/kernel-contract.md @@ -447,7 +447,10 @@ bounds read off the operator (`translate.radial.radial_rate_bounds`): a constant table's exact extremes of `Re D`, and, for a varying medium, a sufficient bound on each interval of the rate's Chebyshev series (`|T_j| ≤ 1`) for a z table or for an expansion's scalar weights against its terms. The refusal names the radial operator's -attenuation and the admissible step. The modal kinds are not checked. +attenuation and the admissible step. The modal kinds are not checked: the default +`LossPolicy` clamps a mode's power attenuation on `D` to `2 × 708 / L`, so no step within +the run forms an exponent past 708; gain, a raised `max_alpha` and the absorber's rates, +added after the clamp, are outside that bound. The imaginary part of Φ is also the unitary phase for the noise field. Modal Φ includes the real endpoint term `-1/2 [ln w(z) - ln w(0)]`. It vanishes identically diff --git a/docs/guide/numerics/boundaries.md b/docs/guide/numerics/boundaries.md index 570b8a773..85a70d9fe 100644 --- a/docs/guide/numerics/boundaries.md +++ b/docs/guide/numerics/boundaries.md @@ -330,7 +330,8 @@ and nothing else. So: long drift through a narrow collar is better served by a wider window. A linear packet walking through this collar at the stated drift leaves `exp(−2Γ)` of its -power, a millionth, to within the one sample the sizing's integral leaves out of the period. +power, a millionth: the sizing sums the rate over the periodic axis, which is exact for this +taper, and the measured exponent agrees with it to a few parts in `10⁹`. The sizing assumes a steady drift over the run, as `warn_walkoff` does. Walk-off that starts late, after fission or as a soliton decelerates, crosses faster than the mean and is absorbed less, so the [monitor](#watching-the-windows-edges) stays the measurement. The default @@ -582,12 +583,27 @@ the spatial routes' default, and about half off the in-stage term, on the same g light can lose most of its band while the run's energy moves by `10⁻⁵`, and the band's loss is what says the window is too short. The rate is diagonal in `ω`, so each bin loses `α(ω) S(ω, z)`, with `S` the `beam_spectrum`, integrated in `z`. That statistic is taken at saved -planes and not [between them](../results/statistics.md#between-the-planes), so the integral runs -over the planes, with an error of about `(α Δz)²/12` of a bin's loss by the trapezoid. Dividing -by what the band kept and lost together counts light the run generated part-way as well as -light it launched. `absorbed_power` has the user's share as its `custom` channel, but over the -whole beam: only this product reads it per band. A lossy medium's own absorption is a further -channel: fused silica absorbs at these wavelengths, and what it takes appears in neither count. +planes and not [between them](../results/statistics.md#between-the-planes), so the integral has +to be read from the planes. + +- **Not by the trapezoid.** Its error, about `(a Δz)²/12` of a bin's loss for a bin falling at + the rate `a`, is small only while `a Δz` is. The bins this rate exists for have `a Δz` large: + light already in one is gone within a fraction of the interval, and the trapezoid charges it + at the plane's value over the whole of it. On a 5 µm slab of fused silica with a 1 fs launch + at 260 nm it overcounted a band's loss by 1.5 to 2 times, and converged slowly as planes + were added. +- **Between two planes.** A bin obeys `dS/dz = G − a S`, with `a` every rate that drains it, + the absorber's and the medium's, and `G` whatever else feeds it. Taking `G` constant over the + interval, the two planes fix it: + `G Δz = x (S₁ − S₀ e⁻ˣ)/(1 − e⁻ˣ)` with `x = a Δz`, and the drains removed `G Δz − (S₁ − S₀)`. + That is exact for a bin that only decays and for one fed at a steady rate, and each rate takes + its share `α/a` of it. + +Dividing by what the band kept and lost together counts light the run generated part-way as +well as light it launched. `absorbed_power` has the user's share as its `custom` channel, but +over the whole beam: only this product reads it per band. A lossy medium's own absorption is a +further channel: fused silica absorbs at these wavelengths, so its rate is one of the drains, +and what it takes is its own share, in neither the walk-off's count nor the collar's. ```python planes = np.linspace(length / 32, length, 32) @@ -611,14 +627,25 @@ output = Simulation.from_free( # (n_z, n_omega), J/(rad/s) per metre of y spectrum = np.asarray(output["stats/beam_spectrum"]) z = np.asarray(output["z"]) -lost = np.trapezoid(walkoff_rate * spectrum, z, axis=0) -collar_lost = np.trapezoid( - (boundary.spectral_rate - walkoff_rate) * spectrum, z, axis=0 -) +glass = np.zeros(grid.n_omega) # the medium's own power loss, 1/m +glass[active] = 2 * grid.omega[active] * np.imag(medium.index(grid.omega[active])) / C +drain = boundary.spectral_rate + glass # every rate that drains a bin +x = drain * np.diff(z)[:, np.newaxis] +# x / (1 - exp(-x)), which is 1 where nothing drains +weight = np.divide(x, -np.expm1(-x), out=np.ones_like(x), where=x > 0) +fed = weight * (spectrum[1:] - np.exp(-x) * spectrum[:-1]) # G dz between the planes +drained = np.sum(fed - np.diff(spectrum, axis=0), axis=0) +per_rate = np.divide(drained, drain, out=np.zeros_like(drained), where=drain > 0) +lost = walkoff_rate * per_rate +collar_lost = (boundary.spectral_rate - walkoff_rate) * per_rate +glass_lost = glass * per_rate wavelength = 2 * np.pi * C / np.where(grid.omega > 0, grid.omega, np.inf) band = (wavelength >= 150e-9) & (wavelength <= 190e-9) share = lost[band].sum() / ( - lost[band].sum() + collar_lost[band].sum() + spectrum[-1, band].sum() + lost[band].sum() + + collar_lost[band].sum() + + glass_lost[band].sum() + + spectrum[-1, band].sum() ) print(f"the walk-off rate took {share:.2f} of the 150-190 nm band") ``` diff --git a/docs/guide/numerics/linear-operator.md b/docs/guide/numerics/linear-operator.md index a1e645c18..ab101292a 100644 --- a/docs/guide/numerics/linear-operator.md +++ b/docs/guide/numerics/linear-operator.md @@ -216,6 +216,13 @@ is never clamped to zero. A non-finite index raises rather than passing through. itself, before $\Phi$ is integrated — clamping the accumulated exponent instead would make the answer depend on the step. +The same clamp keeps a modal run inside the exponent budget that a separable or radial +run is checked against ([free-space model](../free-space/model.md#exact-longitudinal-operator)): no step within the +fibre's length forms a factor past $\mathrm{e}^{708}$, below the overflow at +$\ln(\mathrm{DBL\_MAX}) \approx 709.78$, so a replayed step cannot overflow and +translation does not check it. Gain, a `max_alpha` raised above the default and the +absorber's rates, which are added to `D` after the clamp, are outside that bound. + ## The operator of a dispersion-specified fibre ```{eval-rst} diff --git a/src/solprop/boundaries.py b/src/solprop/boundaries.py index d30f08430..5f7f62a2c 100644 --- a/src/solprop/boundaries.py +++ b/src/solprop/boundaries.py @@ -579,8 +579,10 @@ def _walkoff_collar_length(grid: TimeGrid, walkoff: float) -> float: \ell_c = \frac{L}{\Delta\tau\,\Gamma}\,\frac{A}{2}\int (1 - W)\,dt. The integral is taken as :func:`warn_walkoff` takes it, so that the two agree on any - grid. For two edges of width ``w`` the Planck taper's symmetry makes it ``w``, - giving ``ell_c = 6.3 w L / (Gamma dtau)``. + grid: the rectangle sum over the samples, exact on the periodic axis, where a + trapezoid would drop the one sample that closes the period. For two edges of width + ``w`` the Planck taper's symmetry makes it ``w``, giving + ``ell_c = 6.3 w L / (Gamma dtau)``. - ``walkoff``: the drift over ``grid.length``, s. - Returns: the collar length in metres. @@ -592,7 +594,7 @@ def _walkoff_collar_length(grid: TimeGrid, walkoff: float) -> float: ) raise GridError(msg) unit_rate = temporal_collar_rate(grid.t, grid.t_collar, 1.0) - integral = float(np.trapezoid(unit_rate, dx=grid.dt)) + integral = float(np.sum(unit_rate)) * grid.dt length = ( 0.5 * integral * grid.length / (MIN_COLLAR_ATTENUATION * walkoff) ) * _WALKOFF_ROUNDING_MARGIN @@ -689,7 +691,7 @@ def absorber( >>> bool(np.array_equal(walking.spectral_rate, boundary.spectral_rate)) True >>> round(walking.in_stage_length / boundary.in_stage_length, 3) - 1.023 + 1.025 """ if walkoff is not None and not temporal: msg = ( @@ -1372,9 +1374,9 @@ def warn_walkoff( if group_delay_span <= 0.0: return rate = temporal_collar_rate(grid.t, grid.t_collar, length) - crossing = ( - 0.5 * float(np.trapezoid(rate, dx=grid.dt)) * grid.length / group_delay_span - ) + # The rectangle sum, not a trapezoid: the axis is periodic, and the sum is exact + # for the taper the packet crosses. + crossing = 0.5 * float(np.sum(rate)) * grid.dt * grid.length / group_delay_span if crossing >= MIN_COLLAR_ATTENUATION: return warnings.warn( diff --git a/src/solprop/translate/problem.py b/src/solprop/translate/problem.py index d4c578d60..77d3f79c7 100644 --- a/src/solprop/translate/problem.py +++ b/src/solprop/translate/problem.py @@ -1430,8 +1430,16 @@ def _check_exponent_budget(self, stepper: StepperSettings) -> None: the budget only at its floor ``dz_min``; above it, its controller rejects a non-finite trial and shrinks. The radial kind carries a dense operator and forms the same inverse factor from it, so it is bounded by - :func:`solprop.translate.radial.radial_rate_bounds` instead (SEP 22, Q16); - the modal kinds are not checked. + :func:`solprop.translate.radial.radial_rate_bounds` instead (SEP 22, Q16). + + The modal kinds are not checked because they need no check: a mode's loss + reaches ``D`` through :func:`solprop.linear.dispersion_operator`, whose default + :class:`solprop.linear.LossPolicy` clamps the power attenuation to + ``2 * MAX_ATTENUATION_EXPONENT / length``, so no step within the run's length + forms an amplitude exponent past 708, below ``log(DBL_MAX)``. Outside that + reason, and so unchecked here: gain, which the policy lets through; a + ``max_alpha`` raised above the default limit; and the absorber's rates, which + are added to ``D`` after the clamp. """ operator = self.linear.separable geometry = self.transform.geometry diff --git a/tests/test_boundaries.py b/tests/test_boundaries.py index 20b1c852a..136f6c89c 100644 --- a/tests/test_boundaries.py +++ b/tests/test_boundaries.py @@ -454,9 +454,12 @@ def test_a_malformed_absorber_is_rejected(field_grid): def crossing_nepers(grid, group_delay_span, boundary_n=DEFAULT_BOUNDARY_N): - """The field attenuation a walking packet accumulates crossing the collar.""" + """The field attenuation a walking packet accumulates crossing the collar. + + The rate summed over the periodic axis, which the packet crosses whole. + """ rate = absorber(grid, boundary_N=boundary_n).temporal_rate - return 0.5 * float(np.trapezoid(rate, dx=grid.dt)) * grid.length / group_delay_span + return 0.5 * float(np.sum(rate)) * grid.dt * grid.length / group_delay_span def test_the_walk_off_warning_fires_when_the_absorber_is_too_weak(field_grid): @@ -468,6 +471,32 @@ def test_the_walk_off_warning_fires_when_the_absorber_is_too_weak(field_grid): warn_walkoff(field_grid, span) +@pytest.mark.parametrize("grid_name", GRIDS) +def test_the_walk_off_warning_reads_the_periodic_sum_of_the_collar(grid_name, request): + """The threshold is the closed form ``6.3 w L / (Gamma ell)``, not one sample short. + + Over both edges the rate integrates to ``12.6 w / ell`` exactly on the periodic + axis; a trapezoid over the samples would read ``w - dt`` and fire at a span + ``dt / w`` short of the threshold. + """ + grid = request.getfixturevalue(grid_name) + width = edge_width(grid) + reference = grid.length / DEFAULT_BOUNDARY_N + threshold = ( + 0.5 + * TEMPORAL_ALPHA_LENGTH + * width + * grid.length + / (reference * MIN_COLLAR_ATTENUATION) + ) + assert grid.dt / width > 1e-4 + with warnings.catch_warnings(): + warnings.simplefilter("error") + warn_walkoff(grid, threshold * (1.0 - 1e-6)) + with pytest.warns(SolpropWarning, match="wrapping around"): + warn_walkoff(grid, threshold * (1.0 + 1e-6)) + + def test_the_walk_off_warning_is_silent_when_the_absorber_is_strong_enough(field_grid): span = 100e-15 assert crossing_nepers(field_grid, span) > MIN_COLLAR_ATTENUATION @@ -554,20 +583,15 @@ def edge_width(grid): def test_a_walkoff_sets_the_length_the_derivation_gives(field_grid): """``ell_c = 6.3 w L / (Gamma dtau)``, on the integral ``warn_walkoff`` reads. - The integral of ``1 - W`` over both edges is ``w`` by the Planck taper's symmetry; - the trapezoid on the axis, which stops one sample short of closing the period, - reads ``w - dt``. So the bound is that one sample, and the sampled form is exact. + The integral of ``1 - W`` over both edges is ``w`` by the Planck taper's symmetry, + and the rectangle sum over the periodic axis reads it exactly, so the bound is the + rounding margin of the sizing. """ drift = 3e-12 width = edge_width(field_grid) collar = absorber(field_grid, walkoff=drift) closed_form = 6.3 * width * field_grid.length / (MIN_COLLAR_ATTENUATION * drift) - assert collar.in_stage_length == pytest.approx( - closed_form, rel=2.0 * field_grid.dt / width - ) - assert collar.in_stage_length == pytest.approx( - closed_form * (width - field_grid.dt) / width, rel=1e-6 - ) + assert collar.in_stage_length == pytest.approx(closed_form, rel=1e-8) assert float(collar.temporal_rate.max()) == pytest.approx( TEMPORAL_ALPHA_LENGTH / collar.in_stage_length, rel=1e-3 ) @@ -815,11 +839,12 @@ def test_a_walkoff_collar_removes_light_of_the_stated_walk_off(backend, stretch) """A packet of the stated drift leaves a millionth of its power, and no more. The packet crosses both edges rigidly, so its power falls by exactly the integral - of the rate over the ring divided by the speed: ``exp(-2 Gamma)`` of the sizing, - corrected by the one sample the sizing's trapezoid leaves out of the periodic sum. + of the rate over the ring divided by the speed: ``exp(-2 Gamma)``, which the sizing + reads from the same periodic sum. A drift of half the collar's or twice its stated value shows the exponent is linear in the speed, and that the light has no other way out. The bound on the - exponent is the controller's ``rtol``; measured 6e-9 on both backends. + exponent is the controller's ``rtol``; measured 2.4e-9 at the stated drift on + ``numpy_ref``. """ if backend == "native": pytest.importorskip("solprop_native") @@ -830,7 +855,7 @@ def test_a_walkoff_collar_removes_light_of_the_stated_walk_off(backend, stretch) collar = absorber(grid, walkoff=drift / stretch) closed_period = np.sum(temporal_collar_rate(grid.t, grid.t_collar, 1.0)) * grid.dt expected = closed_period * grid.length / (collar.in_stage_length * drift) - assert expected == pytest.approx(2.0 * MIN_COLLAR_ATTENUATION / stretch, rel=0.02) + assert expected == pytest.approx(2.0 * MIN_COLLAR_ATTENUATION / stretch, rel=1e-8) power = transmitted_power(grid, collar, drift, backend=backend) assert -math.log(power) == pytest.approx(expected, rel=1e-7) @@ -846,8 +871,7 @@ def test_without_a_temporal_collar_the_same_packet_loses_nothing(): def test_the_split_scheme_removes_the_same_light(): """The split's factors telescope, within its crossing limit, to the same exponent. - The split's first-order error in the step and the sizing's left-out sample, 0.7 % - here, set the bound; measured 0.5 %. + The split's first-order error in the step sets the bound; measured 0.23 %. """ grid = EnvelopeGrid.from_limits(1.0, 800e-9, (600e-9, 1200e-9), 2e-12, thg=False) drift = 0.8 * grid.n_t * grid.dt diff --git a/tests/test_spectral_walkoff_recipe.py b/tests/test_spectral_walkoff_recipe.py index e1c3936f7..d33e5beca 100644 --- a/tests/test_spectral_walkoff_recipe.py +++ b/tests/test_spectral_walkoff_recipe.py @@ -294,6 +294,56 @@ def test_the_per_band_loss_read_from_planes_is_the_loss(backend): ) +def interval_drained(spectrum, z, drain): + """What every drain took from each bin, read exactly between the planes. + + The guide's formula: with ``dS/dz = G - a S`` and ``G`` constant over an interval, + ``G dz = x (S1 - S0 exp(-x)) / (1 - exp(-x))``, ``x = a dz``, and the drains took + ``G dz - (S1 - S0)``. Returns that sum over the intervals, divided by ``a``, so a + rate times it is the rate's share. + """ + x = drain * np.diff(z)[:, np.newaxis] + weight = np.divide(x, -np.expm1(-x), out=np.ones_like(x), where=x > 0) + fed = weight * (spectrum[1:] - np.exp(-x) * spectrum[:-1]) + drained = np.sum(fed - np.diff(spectrum, axis=0), axis=0) + return np.divide(drained, drain, out=np.zeros_like(drained), where=drain > 0) + + +def test_the_guide_s_interval_formula_is_exact_on_few_planes(): + """The bins the rate is for fall far between planes; the formula still reads them. + + On four planes ``a dz`` reaches 2.7 nepers at the launch's bins, where the + trapezoid overcounts. The interval formula is exact for a bin that only decays, so + the controller's ``rtol`` and rounding are the bound, and four planes and + thirty-two agree to it. + """ + rate = recipe_rate(math.log(1 / SURVIVAL) / LENGTH) + glass = np.zeros(GRID.n_omega) + active = GRID.active_mask + glass[active] = ( + 2 * GRID.omega[active] * np.imag(MEDIUM.index(GRID.omega[active])) / C + ) + whole = rate + glass + field = launch([(carrier_at(ONSET + 0.5 * PROTECTED), 1e-9, 0.0)]) + readings = [] + for planes in (4, 32): + output = run(spectral_only(rate), field, planes=planes) + spectrum = np.asarray(output["stats/beam_spectrum"]) + z = np.asarray(output["z"]) + readings.append(rate * interval_drained(spectrum, z, whole)) + few, many = readings + weak = (spectrum[0] > 1e-3 * spectrum[0].max()) & (rate > 0.0) + assert np.count_nonzero(weak) > 3 + assert float(np.max(whole[weak])) * LENGTH / 4 > 2.0 + exact = ( + spectrum[0, weak] + * (rate[weak] / whole[weak]) + * -np.expm1(-whole[weak] * LENGTH) + ) + np.testing.assert_allclose(few[weak], exact, rtol=1e-6) + np.testing.assert_allclose(few[weak], many[weak], rtol=1e-6) + + def test_a_band_can_lose_most_of_itself_while_the_run_s_energy_barely_moves(): """Why the check is per band: the total reads ``1e-5`` while the band reads tens of per cent.""" diff --git a/tests/translate/test_operator_settings.py b/tests/translate/test_operator_settings.py index 9806e6778..335dc243e 100644 --- a/tests/translate/test_operator_settings.py +++ b/tests/translate/test_operator_settings.py @@ -25,7 +25,8 @@ TranslationError, ) from solprop.grid import EnvelopeGrid, LineGrid, RadialGrid -from solprop.linear import Frame +from solprop.linear import MAX_ATTENUATION_EXPONENT, Frame +from solprop.nonlinear.gnlse import SpecifiedMode from solprop.nonlinear.responses import KerrResponse from solprop.output.serialisation import serialise from solprop.profiles import Linear @@ -504,3 +505,61 @@ def test_the_radial_bounds_cover_a_varying_medium(): # And close: the series converge fast on a smooth density, so the bound # sits within a hundredth of the sampled extreme. assert bounds.loss <= 1.01 * float(np.max(-sampled)) + + +# ------------------------------------------------------------ the modal kinds + + +def lossy_fibre(alpha: float, integrator) -> Simulation: + """A 4 cm fibre of power loss ``alpha``, replayed in four steps of 1 cm.""" + length = 0.04 + grid = EnvelopeGrid.from_limits(length, WAVELENGTH, (600e-9, 1.2e-6), 2e-12) + mode = SpecifiedMode.from_taylor( + [0.0, 0.0, -2e-26], omega_ref=grid.omega0, loss=alpha, aeff=8e-11 + ) + pulse = GaussianPulse(wavelength=WAVELENGTH, fwhm=100e-15, energy=1e-9) + step = length / 4 + return Simulation( + grid, + mode, + [], + pulse, + settings=SimulationSettings( + controller=PIController(rtol=1e-6, atol=1e-14), + save_planes=[length], + dz_max=step, + dz_init=step, + schedule=FixedSchedule.uniform(0.0, length, step, save_planes=[length]), + integrator=integrator, + backend="numpy_ref", + absorber=bare_absorber(grid), + boundary_check=None, + prefactor_index=1.0, + ), + noise=None, + ) + + +@pytest.mark.parametrize( + "integrator", [RungeKutta(), Multistep()], ids=["runge_kutta", "multistep"] +) +def test_the_loss_clamp_keeps_a_modal_replay_inside_the_budget(integrator): + """The modal kinds go unchecked because ``LossPolicy`` bounds their exponent. + + The loss asks for ``alpha h / 2 = 5000`` nepers a step, past ``log(DBL_MAX)``; + the clamp holds the rate to ``2 * 708 / L``, a quarter of 708 a step of + ``L / 4``, so translation neither refuses nor warns and the field stays finite. + """ + alpha = 1e6 + with warnings.catch_warnings(): + warnings.simplefilter("error", ExponentBudgetWarning) + simulation = lossy_fibre(alpha, integrator) + problem = simulation.translate() + step = simulation.grid.length / 4 + assert alpha * step / 2 > LOG_DBL_MAX + rate = problem.linear.phi_constant + assert rate is not None + exponent = float(np.max(-np.real(rate))) * step + assert exponent == pytest.approx(MAX_ATTENUATION_EXPONENT / 4, rel=1e-12) + output = simulation.run(progress="none") + assert np.all(np.isfinite(output["field"]))