fdtd: waveguide modes, mode monitor and source, and Gaussian beams (stage 7)

docs/_toc.yml

@@ -113,6 +113,8 @@ subtrees:
           - file: guide/fdtd/general-fdtd
           - file: guide/fdtd/model
           - file: guide/fdtd/time-domain
+          - file: guide/fdtd/waveguide-modes
+          - file: guide/fdtd/gaussian-beam
   - file: reference/index
     title: API reference
     subtrees:
@@ -252,6 +254,10 @@ subtrees:
               - file: validation/fdtd/fd-public-numerics-refusals
               - file: validation/fdtd/fd-subpixel-occupancy
               - file: validation/fdtd/fd-two-debye-film
+              - file: validation/fdtd/fd-waveguide-modes
+              - file: validation/fdtd/fd-waveguide-through-pml
+              - file: validation/fdtd/fd-mode-source
+              - file: validation/fdtd/fd-gaussian-beam
               - file: validation/fdtd/fdtd-adjoint-validation
               - file: validation/fdtd/fdtd-bragg-cavity
               - file: validation/fdtd/fdtd-circle-gradients

docs/design/unified-api.md

@@ -368,7 +368,7 @@ claim them without changing the adapter and geometry accounting.
 
 The standard FDTD problem can request time records, fields, flux, diffraction,
 far fields, and final fields. Only implemented observables are public stable
-exports. GaussianBeam, mode sources/monitors, decay termination, GPU support,
+exports. Decay termination, GPU support,
 and nonuniform grids must not appear as supported merely because constructors
 or keywords have been reserved.
 
@@ -1194,7 +1194,10 @@ not extend that scope.
 | `rcwa.operators` and `rcwa.experimental` (§5) | implemented, as submodules imported explicitly | `src/wavebit/rcwa/operators.py`, `src/wavebit/rcwa/experimental/` |
 | Detuning as explicit, recorded policy (§8) | implemented: `rcwa.Numerics(detuning=None)` is the default and solves the requested point | `src/wavebit/rcwa/numerics.py` |
 | `fdtd.Problem`, `Scene`, `Domain`, `Boundaries`, `Structure`, `PlaneWave`, `GaussianPulse`, `Outputs`, `Record`, `Numerics`, `UniformGrid`, `Duration`, `PML`, with the PML outside the physical domain (§6) | implemented; CPU and double precision only (§13) | `src/wavebit/fdtd/` |
-| FDTD `GaussianBeam`, mode sources and monitors, decay termination, nonuniform grids, GPU (§6) | not implemented, and not exposed | -- |
+| FDTD decay termination, nonuniform grids, GPU (§6) | not implemented, and not exposed | -- |
+| FDTD discrete waveguide modes (`solve_modes`) and `ModeMonitor` (stage 7) | implemented | [validation](../validation/fdtd/fd-waveguide-modes.md) |
+| FDTD `ModeSource`, a one-way mode launch (stage 7) | implemented | [validation](../validation/fdtd/fd-mode-source.md) |
+| FDTD `GaussianBeam`, a focused beam launched one way (stage 7) | implemented | [validation](../validation/fdtd/fd-gaussian-beam.md) |
 | `fdtd.Numerics.geometry_ranges` (§7) | implemented | `src/wavebit/fdtd/_numerics.py` |
 | Host-only topology decisions rejected under an external transform (§7) | implemented: fixed-cell Cartesian Li rectangle sizes/centres use an array-valued edge-order guard under JIT/vmap. Other scalar region geometry can explicitly select `topology_check='host'` for checked CPU callbacks. The default device policy refuses it under `jax.jit` or `jax.vmap` (`prepared.host_topology`); an eager `jax.grad` tracer's primal is checked as a PyTorch tensor is | `src/wavebit/rcwa/_lifecycle.py` |
 | `outputs`, `validity`, `diagnostics`, `provenance` on every result; `DiffractionOrders` with `order(m, n)`, `FieldSamples`, `SpectralFieldSamples`, `RawDiffractionSpectrum`, `TimeSeries` (§8) | implemented | `src/wavebit/observations/` |

docs/development/source-map.md

@@ -86,6 +86,8 @@ research cases; inspect its scope before launching a long calculation.
 | Coherent sources and absolute power | [test_unified_rcwa_sources.py](https://github.com/fieldalgo/wavebit/blob/main/tests/api/test_unified_rcwa_sources.py) | SI source strength, finite observation regions and independently evaluated fields. |
 | FDTD domain, source/monitor or reuse | [test_lifecycle.py](https://github.com/fieldalgo/wavebit/blob/main/tests/fdtd/test_lifecycle.py) | PML outside the physical domain, unchanged laboratory coordinates, live source/material reuse and resolved duration guards. |
 | RCWA modal or reconstruction algebra | [test_smm_field_state.py](https://github.com/fieldalgo/wavebit/blob/main/tests/rcwa/test_smm_field_state.py), [test_lossy_media.py](https://github.com/fieldalgo/wavebit/blob/main/tests/rcwa/test_lossy_media.py) | Independent interior fields and interface responses rather than a regenerated solver golden. |
+| FDTD waveguide modes | [test_fd_waveguide_modes.py](https://github.com/fieldalgo/wavebit/blob/main/tests/fdtd/test_fd_waveguide_modes.py) | Modes of the discrete Yee system against closed forms, exact slab dispersion and a published corner benchmark, and a solved mode advanced by the time stepper itself. |
+| FDTD mode source | [test_fd_mode_source.py](https://github.com/fieldalgo/wavebit/blob/main/tests/fdtd/test_fd_mode_source.py) | A one-way mode sheet's injection and leakage against their exact grid predictions, absolute power, convergence, a butt joint against mode matching, backends, derivatives and refusals. |
 | FDTD time step or constitutive update | [test_fd_stage3_backends.py](https://github.com/fieldalgo/wavebit/blob/main/tests/fdtd/test_fd_stage3_backends.py), [test_fd_constitutive_absorber.py](https://github.com/fieldalgo/wavebit/blob/main/tests/fdtd/test_fd_constitutive_absorber.py) | Unchanged raw fields under observation-only changes, discrete energy/work identities and backend agreement. |
 | Cross-solver phase, port or diffraction normalization | [test_public_lifecycle_diffraction.py](https://github.com/fieldalgo/wavebit/blob/main/tests/cross_solver/test_public_lifecycle_diffraction.py) | Current public APIs, unequal structures and complex nonzero orders; equal total power alone can hide a phase error. |
 | Recorded results and source identity | [test_run_records.py](https://github.com/fieldalgo/wavebit/blob/main/tests/api/test_run_records.py) | Exact arrays/units, corrupt or incomplete archive refusal, clean input association and loading without the producing backend. |

docs/guide/fdtd/gaussian-beam.md

@@ -0,0 +1,59 @@
+# Gaussian beams
+
+`fdtd.GaussianBeam` launches a focused Gaussian beam from a plane, one way.
+
+```python
+from wavebit import Units, fdtd
+
+beam = fdtd.GaussianBeam('b', plane=fdtd.Plane(z=-4.), direction='+z', waist=1.5,
+                         focus=(0., 0., 0.), polarization='y',
+                         pulse=fdtd.GaussianPulse(wavelength=1., frequency_sigma_hz=2e13))
+```
+
+## What it launches
+
+On the focal plane, the E component along `polarization` (a transverse axis) is
+`exp(-r**2 / waist**2)` about `focus`, and the other transverse component is zero.
+`waist` is the 1/e amplitude radius, as for `rcwa.GaussianBeam`. The beam is the propagating
+part of that plane's angular spectrum. It is built from the grid's own plane waves in the
+homogeneous medium that fills the source plane, so at the pulse's centre wavelength it is an
+exact solution of the grid. `focus` is where the beam would focus if that medium filled
+space. A beam that enters another medium focuses where refraction takes it.
+
+`amplitude` is in square-root watts (per metre of each invariant axis): a continuous wave of
+amplitude `A` carries `A**2` watts. `phase` is in radians, and `delay_seconds` delays the
+pulse. In a two-dimensional problem (one `InvariantAxis`) the beam is a sheet beam, focused
+in the resolved transverse axis only.
+
+The plane is a one-sided total-field/scattered-field sheet, as for `ModeSource` (see
+[Waveguide modes](waveguide-modes.md)). At the centre wavelength it launches exactly the
+beam, one way. Away from it, each plane wave's profile and half-cell phase are the centre's,
+and a little of the beam goes backward. The [validation record](../../validation/fdtd/fd-gaussian-beam.md)
+measures that part against its exact prediction.
+
+## Accuracy on a finite grid
+
+The beam is exact for the grid, not for the continuum. The grid's plane waves travel at a
+phase velocity that depends on their direction, so the phase a beam gathers through its focus
+(its Gouy phase) differs a little from the continuum beam's. In the
+[validation record](../../validation/fdtd/fd-gaussian-beam.md) (vacuum, wavelength 1 µm,
+waist 1.5 µm) it differs by at most 2.2e-3 rad within 3 µm of the focus at a spacing of
+0.05 µm (20 cells per wavelength), and by 5.4e-4 rad at 0.025 µm: the error falls at second
+order. Converge the spacing for any result that depends on the phase along the beam. The
+focal-plane waist, the power and the back-radiation match their predictions far more closely.
+
+## What it requires
+
+- The faces across the plane are `Open` (or `InvariantAxis`). A metal, mirror or periodic
+  face is refused (`source.boundary`): the beam would become the modes of that guide or
+  lattice.
+- The medium on the plane is uniform (`source.medium`), lossless and non-dispersive, and the
+  same along the axis for a cell each way (the window rules of `solve_modes`).
+- The cross-section holds the beam: on the outermost cells of the plane it must be at most
+  1e-3 of its peak (`source.beam_size`). The sheet stops at the window, one cell inside an
+  `Open` face, and launches nothing beyond it.
+- `waist` and `focus` are concrete numbers (`source.traced`): the beam is built on the host.
+  `amplitude`, `phase` and `delay_seconds`, and structures elsewhere, may carry traced values.
+  The pulse wavelength may not (`mode.traced`).
+- Normal incidence only: the beam's axis is the plane's normal. A flux monitor within one
+  cell of the sheet is refused (`monitor.source_sheet`).

docs/guide/fdtd/index.md

@@ -15,3 +15,5 @@ Start with the [current FDTD problem API](general-fdtd.md), [how the solver mode
 - [Current FDTD domains, sources and observations](general-fdtd.md)
 - [How the FDTD solver models a problem](model.md)
 - [A solve in the time domain](time-domain.md)
+- [Waveguide modes of the discrete grid](waveguide-modes.md)
+- [Gaussian beams](gaussian-beam.md)

docs/guide/fdtd/model.md

@@ -207,7 +207,7 @@ The rule has limits, each enforced by `check`:
 | The smallest extent the body bounds must be at least 4 cells; 4 to 8 cells warns | across fewer cells, where the mixed cells fall decides the answer more than the body does |
 | Material: `Material` with a scalar permittivity and permeability; no per-axis values | the mixing rule is applied per component and assumes an isotropic body |
 | No `ConductiveMaterial` or `DispersiveMaterial`, except a flat film spanning the cell | lossy and dispersive mixing is validated only for those film scopes |
-| Must not touch an `Open()` face unless an unturned, transversely uniform `Box` explicitly declares `Structure.continuation` | finite-body continuation is not inferred; see [smoothed substrates](substrate-gradients.md) |
+| Must not touch an `Open()` face unless an unturned, transversely uniform `Box` explicitly declares `Structure.continuation` | finite-body continuation is not inferred; see [smoothed substrates](substrate-gradients.md). A cell-snapped waveguide core may continue instead (see [General-scene limits](#general-scene-limits)) |
 | Must not reach a `Periodic()` face, unless it runs the whole period reaching **more than** one cell past both faces | otherwise part of it would be lost at the seam |
 | Two smoothed surfaces may not pass through the same cell | the mixing rule describes two media and one surface per cell; separate them by a cell or write one body |
 | A smoothed surface may not cut a cell that a cell-snapped face also cuts | the cell would hold three media (see below) |
@@ -329,7 +329,7 @@ type per axis or a `(low, high)` pair.
 
 | Face | Meaning |
 | --- | --- |
-| `Open()` | outgoing: a PML outside the physical domain, matched to the uniform medium that reaches the face. Only uniform media may reach it — the background, a `SnappedSlab` across that axis, or an unturned `Box` that declares `Structure.continuation` and reaches more than one cell past both transverse faces; any other finite body touching it, or a slab or continued `Box` crossing an orthogonal `Open()` face, is refused (see [General-scene limits](#general-scene-limits)). |
+| `Open()` | outgoing: a PML outside the physical domain, matched to the uniform medium that reaches the face. Only uniform media may reach it — the background, a `SnappedSlab` across that axis, an unturned `Box` that declares `Structure.continuation` and reaches more than one cell past both transverse faces, or a `SnappedBox` core that declares continuation through that face; any other finite body touching it, or a slab or continued `Box` crossing an orthogonal `Open()` face, is refused (see [General-scene limits](#general-scene-limits)). |
 | `Periodic()` | both faces wrap with no phase; normal incidence in a periodic cell. |
 | `Bloch(phase=φ)` or `Bloch(k=k)` | both faces wrap with a phase; the run becomes complex (twice the memory and work). `Bloch.from_source(source, units=...)` writes the faces an oblique `PlaneWave` needs. |
 | `PEC()` / `PMC()` | tangential E / tangential H vanishes on the wall. A structure reaching past one is clipped, with `GeometryClippedWarning`. |
@@ -389,7 +389,8 @@ These messages do not relax the return tolerance or the measured support limits.
 
 ## Sources
 
-`Problem.sources` accepts `PlaneWave`, `CurrentDipole` and `TFSF` declarations with
+`Problem.sources` accepts `PlaneWave`, `CurrentDipole`, `TFSF`, `ModeSource` and
+`GaussianBeam` declarations with
 unique, nonempty string names and a `GaussianPulse`. `check` and
 `Prepared.validate` report malformed declarations at `sources[i]`, invalid or
 duplicate names at `sources[i].name`, and invalid pulses at
@@ -422,6 +423,8 @@ one matching the cell boundaries, is unsupported and produces a structured
 | `PlaneWave(name, plane=..., pulse=..., direction="-z", theta=0, phi=0, polarization="s")` | a one-way sheet at the node plane nearest `plane`, launching toward `direction`. `theta` is the polar angle from the propagation axis and `phi` the azimuth from the first transverse axis (x for a z wave); `"s"` has E perpendicular to the plane of incidence, `"p"` in it. The transverse faces must be `Periodic()` (normal incidence) or the `Bloch` faces from `Bloch.from_source`; the transverse wavevector is fixed at the pulse's centre wavelength, so an oblique broadband pulse is not constant-angle illumination. |
 | `CurrentDipole(name, center=..., pulse=..., component="Ey", kind="electric", current_moment=1)` | a point **current** (A·m), or with `kind="magnetic"` its magnetic equivalent (V·m), spread over a smooth stamp (`"bspline"`, the default, or `"tent"`, each with a half-width of two cells) that keeps the moment and its centroid exact wherever the centre sits. It is not a dipole moment in C·m. |
 | `TFSF(name, box=..., pulse=..., direction="-z", polarization="x")` | a plane wave inside a `SnappedBox` only (total field inside, scattered field outside), for a finite scatterer. Axis directions only; the two faces along its axis must be `Open()`. |
+| `ModeSource(name, plane=..., pulse=..., direction="+z", mode=0, amplitude=1)` | a waveguide mode of `solve_modes`, launched one way from a sheet at the node plane nearest `plane`; see [Waveguide modes](waveguide-modes.md). |
+| `GaussianBeam(name, plane=..., pulse=..., waist=..., focus=(0, 0, 0), polarization="x", direction="+z")` | a focused Gaussian beam launched one way from a sheet at the node plane nearest `plane`; see [Gaussian beams](gaussian-beam.md). |
 | `InitialField(E, H, amplitude=1)` | the field the run starts from: E at `t = 0` and H at `t = -dt/2`, given as arrays or as callables evaluated at each component's own Yee points. |
 
 `PlaneWave.plane` must be an `fdtd.Plane`, for example `fdtd.Plane(z=.5)`;
@@ -709,18 +712,22 @@ coating of 0.105 µm instead, the same scene would warn `geometry.cell_snapped`
   other GPU backends and single precision are refused. The grid is uniform and cubic.
 - Per-cell sampled grids use `fdtd.native`; this general-scene formulation has no raster body.
   Constructive solid geometry (union, intersection, difference) and tilted extrusion axes are unsupported.
-- No Gaussian beams, mode sources or mode monitors.
+- Mode sources, mode monitors and the mode solver are described in
+  [Waveguide modes](waveguide-modes.md), and Gaussian beams in [Gaussian beams](gaussian-beam.md).
 - No smoothed lossy or dispersive bodies outside flat films.
-- **A body narrower than the cross-section cannot continue into an absorbing layer.** Only
-  uniform media continue into an `Open()` face's PML: the background, a `SnappedSlab` across that
-  axis, and an unturned `Box` that declares `Structure.continuation` and reaches more than one cell
-  past both transverse faces (a substrate; see [smoothed substrates](substrate-gradients.md)). A
-  finite body that touches an `Open()` face is refused — a `Cylinder`, `Prism` or `Sphere`, or a
-  `Box` or `SnappedBox` core — and so is a smoothed body that ends so close to the face that its
-  smoothed edge reaches the layer. A waveguide whose core, round or rectangular, continues into
-  the absorber is therefore not supported yet, and no date is set. Ending the core before the
-  absorber, with its smoothed edge clear of the layer, avoids the refusal, but the end face then
-  reflects, so the scene solves a different problem.
+- **Only a cell-snapped core continues into an absorbing layer.** Into an `Open()` face's PML
+  continue the background, a `SnappedSlab` across that axis (or, declaring
+  `Structure(continuation=...)`, through the faces of the axes it runs along), an unturned `Box`
+  that declares `Structure.continuation` and reaches more than one cell past both transverse
+  faces (a substrate; see [smoothed substrates](substrate-gradients.md)), and a `SnappedBox`
+  core that declares `Structure(continuation=('z-', 'z+'))` (any one axis) and reaches those
+  physical faces. The core is extended through the padding with its cross-section unchanged, so
+  every layer of the PML sees the same medium; its measured guided-mode reflection is in the
+  [validation record](../../validation/fdtd/fd-waveguide-through-pml.md). A smoothed core
+  narrower than the cross-section — a `Box`, `Cylinder`, `Prism` or `Sphere` — is still refused
+  at an `Open()` face, and so is a smoothed body that ends so close to the face that its
+  smoothed edge reaches the layer, or a continued core that touches an `Open()` face of another
+  axis.
 
 Custom static geometry implements `fdtd.Shape`; prepared reuse checks its declared instance state.
 Its cell support must stay clear of Open padding. See the [migration guide](../../migration/index.md)