hz-b · simonevadi · Aug 17, 2026 · Aug 15, 2026 · Aug 15, 2026 · Aug 15, 2026
diff --git a/.gitignore b/.gitignore
@@ -110,3 +110,7 @@ __pycache__/
 *.pyc
 .pytes*t_cache/
 .ruff_cache
+
+# Per-solver checkpoint directories written by run_rcwa.py / run_neviere.py.
+# The unsuffixed checkpoints/ baseline stays tracked; these are scratch.
+validation/*/results/checkpoints_*/
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,6 +1,50 @@
 # Changelog
 
 ## Unreleased
+
+- Every example that solves now accepts `--solver rcwa|neviere`, defaulting to `rcwa`. Solver-dependent example outputs are suffixed (`*_rcwa.*`, `*_neviere.*`) so the two runs sit side by side; geometry artifacts such as `*_profile.png` stay unsuffixed because they do not depend on the solver.
+- Added three examples covering the cases where the two solvers genuinely differ, rather than duplicating existing examples whose curves would be indistinguishable: `deep_grating_limits` (the modal solver stops at 8.4 wavelengths of depth where the differential method reaches 167), `continuous_vs_staircase` (continuous z-sampling is bit-identical across every `z_resolution_nm`, while a staircase run carries 2.7e-3 of discretization error at 2 nm), and `solver_runtime` (2.4x to 3.0x faster at production resolution, with the solvers still agreeing to ~1e-11).
+- Fixed `multilayer_theta_search` and `blazed_multilayer_memory_comparison` failing with `BrokenProcessPool` on macOS. Both use `max_workers`, and a spawned worker re-imports the script by path; without a `__main__` guard each worker re-ran the whole example and recursively spawned more.
+
+- Fixed `run_multilayer_theta_search_sweep` reporting `solver="rcwa"` on every case regardless of the solver requested. The sweep builds its own `CaseExecutionResult` and copied only the theta-search diagnostics from the underlying result, so `solver` and `solver_options` fell back to their dataclass defaults. The computation was correct throughout -- an rcwa and a neviere sweep differ by ~5e-14 -- but a saved sweep recorded the wrong provenance.
+- Extended the both-solver test coverage to every entry point the examples use: `monochromator_cases`, `energy_angle_cases`, `run_parameter_study`, `run_multilayer_theta_search_sweep`, `assemble_custom_stack`/`LayerSpec`, `AFMGrating`, and the `write_all_orders_csv`/`plot_order_subset`/`efficiency_for_order` output helpers. Swapping `--solver` on any example now exercises a tested path.
+
+- Polarization arguments now accept `TE` and `TM` alongside `s` and `p`, case-insensitively, through a single shared `grax.normalize_polarization`. The library already spoke TE/TM internally (`res0`, `_solve_te_stack`) and in its theory docs while the public API took only s/p; these are the same two states. Values canonicalize to `s`/`p`, so results, CSVs and checkpoints are unchanged. Note the equivalence holds in classical mounting, which is the only mounting this solver supports.
+- The multilayer theta-search workflow now takes `polarization`. It previously had no such argument anywhere, so every theta search ran the `run_simulation` default of `s` -- awkward given all four validation cases run `p`. `run_multilayer_theta_search`, `run_multilayer_theta_search_sweep`, `multilayer_theta_search_cases` and the web form all accept it, and it reaches all three stages of the search. Default stays `s`.
+- `run_parameter_study` now validates `polarization` at the entry point instead of failing several frames deeper inside the simulation wrapper.
+
+### Breaking changes
+
+Renames only — no behaviour changed. Migration:
+
+| old | new |
+| --- | --- |
+| `RCWASimulation` | `GratingSimulation` |
+| `BatchSimulationRunner(default_diffraction_order=)` | `diffraction_order=` |
+| `BatchSimulationRunner(default_fourier_orders=)` | `fourier_orders=` |
+| `BatchSimulationRunner(default_polarization=)` | `polarization=` |
+| `BatchSimulationRunner(default_solver=)` | `solver=` |
+| `neviere_options=` (all entrypoints) | `solver_options=` |
+| `run_simulation(min_efficiency=)` | `min_reflected_efficiency=` |
+
+- Renamed `grax.simulation.RCWASimulation` to `GratingSimulation`. The class drives whichever solver `solver=` selects, so the old name described only one of the two things it can do.
+- Dropped the `default_` prefix from the `BatchSimulationRunner` settings. The prefix was accurate — each of these is the value a case inherits when it omits its own key — but it made `default_fourier_orders: int = 25` read as the default of a default, and the runner's arguments now mirror `run_simulation`'s exactly. The per-case override relationship is documented on each argument instead.
+- Renamed the `neviere_options` argument to `solver_options` everywhere, matching the `solver_options` field already present on results. One name in both directions, and the generic runner no longer names a specific solver in its API. The `NeviereOptions` class keeps its name.
+- Renamed `min_efficiency` to `min_reflected_efficiency` on `run_simulation`, `GratingSimulation` and `run_multilayer_theta_search`. It was half of a min/max pair whose other half was already `max_reflected_efficiency`, and the batch runner had to translate between the two spellings.
+
+- Added a second electromagnetic solver, the Nevière differential method, selectable with `grax.run_simulation(..., solver="neviere")` and `BatchSimulationRunner(default_solver=...)` (per-case `"solver"` also works). It expands the fields in the same truncated Fourier basis as RCWA and applies the same Li/fast-Fourier-factorization rules for TM, but integrates the coupled first-order system in `z` with fourth-order Runge-Kutta instead of eigen-decomposing each layer. Every default stays `"rcwa"`, so existing code, checkpoints and artifacts are unaffected. References: Nevière, Vincent & Petit, *Nouv. Rev. Optique* **5**, 65 (1974); Nevière, *JOSA A* **11**, 1835 (1994); Nevière & Popov, *Light Propagation in Periodic Media* (CRC, 2003).
+- `grax.NeviereOptions` controls the differential method's integration. Step and sub-block sizes are given in optical phase rather than nanometers, so one setting behaves consistently across photon energies, grazing angles and truncation orders. The Fourier truncation order remains the existing `fourier_orders` argument.
+- `NeviereOptions(z_sampling="continuous")` drops the staircase approximation both solvers otherwise share: the permittivity is re-expanded from the true grating profile every `sample_phase` of optical depth, so the result no longer depends on `z_resolution_nm`.
+- The differential method is numerically more robust on deep gratings. The modal solver evaluates `q / sinh(q d)` across a whole layer and overflows above roughly seven wavelengths of depth for a high-contrast lamellar grating; the differential method caps the optical thickness of anything it forms explicitly and still conserves energy to `1e-9` at 167 wavelengths.
+- `SingleSimulationResult` and `CaseExecutionResult` gained a `solver` field, round-tripped through checkpoints so a resumed sweep keeps that provenance.
+- Split the 1D solver into a `grax.solvers` package: `solvers/common.py` holds the shared types, `res0`/`res1`, the Fourier machinery, the layer field operators, the interface cascade and the efficiency extraction; `solvers/rcwa.py` holds the modal layer solve and `res2`; `solvers/neviere.py` holds the differential method. `grax.rcwa_1d` remains as a re-export shim, and the RCWA numerics are unchanged (verified bit-identical across laminar, blazed, multilayer and sinusoidal cases in both polarizations).
+- The batch progress bar now names the solver actually running (`neviere batch`) instead of always printing `RCWA batch`, and the solver-agnostic docstrings on `SingleSimulationResult`, `BatchSimulationRunner`, `gratings.py` and `run_parameter_study` no longer describe themselves as RCWA-specific.
+- Solver selection now reaches every workflow. `run_multilayer_theta_search`, `run_multilayer_theta_search_sweep`, `run_parameter_study`, the `grax_opt` measurement fits and the web UI run form all accept a solver, alongside the existing `run_simulation` and `BatchSimulationRunner` support. Every default stays `"rcwa"`.
+- Fixed the multilayer theta-search workflow silently ignoring `solver=`. `BatchSimulationRunner(default_solver="neviere")` computed those cases with RCWA, with no error or warning, because the workflow's payload carried `backend` but not `solver`. Both runner-settings mappings now come from one place so a new setting cannot reach one execution path and miss the other.
+- `run_parameter_study` gained `backend`, which it previously hardcoded to `"numba"` internally.
+- Results now record `solver_options` alongside `solver`, so a checkpointed differential-method run pins the integration settings that produced it and not just the solver name.
+- Reorganized each validation case into a `grating_definition.py` holding the grating and sweep grid, plus a `run_rcwa.py` and a `run_neviere.py` that import it. Both solvers therefore see identical geometry, energy grid and truncation by construction. Each case's comparison script now overlays both solver curves alongside the external reference codes.
+- Fixed the validation sweep scripts crashing with `BrokenProcessPool` on macOS: the batch runner spawns workers there, and without a `__main__` guard each spawned worker re-ran the whole sweep and recursively spawned more workers.
 - `random-interface` roughness is now a correlated Gaussian random field (Gaussian autocorrelation) instead of per-sample white noise. `grax.RoughnessSpec` gains `correlation_length_nm`: the lateral autocorrelation length in nanometers, defaulting to one tenth of the grating period (`0.0` reproduces the previous white-noise interface). This produces physically smooth interfaces; long correlation lengths (much larger than the grating period) wash out geometrically and are better modelled by the `debye-waller` kind.
 - Added per-layer roughness: `LayerSpec` accepts `roughness_sigma_nm`, and `SingleLayerStack`/`MultilayerStack`/`CustomStack` expose per-layer/per-material roughness arguments (plus a `substrate_roughness_sigma_nm` for the substrate boundary). Each value sets the roughness of that interface, falling back to the grating-level `RoughnessSpec.sigma_nm` when unset. `random-interface` perturbs each interface with its own sigma; per-layer `debye-waller` sigmas combine in quadrature into an effective damping.
 - Web UI: the grating form now has a roughness sigma field for the substrate, each coating layer/material, and the top cap (persisted with the grating, schema v2). The run form has a roughness-kind dropdown (None / Debye-Waller / Random interface) that selects the model applied at run time.

diff --git a/README.md b/README.md
@@ -1,15 +1,24 @@
-# grax
+# GraxPy
 
-`grax` is an independent Python package for diffraction-grating simulations in
-X-ray optics, inspired by RETICOLO v9 and extended with higher-level workflows
-for practical studies.
+`graxpy` is an independent Python package for one-dimensional X-ray
+diffraction-grating simulations. Its public Python import is `grax`.
+
+`grax` provides two selectable electromagnetic solver paths:
+
+- modal rigorous coupled-wave analysis (RCWA), inspired by RETICOLO v9;
+- the Nevière differential method.
+
+Both solvers support the same grating, material, and workflow APIs. They share
+the Fourier/discretization infrastructure and differ in how they propagate the
+fields through a layer.
 
 ## Documentation
 
-Full user and API documentation is published online.
+Full user and API documentation is available at
+[graxpy.readthedocs.io](https://graxpy.readthedocs.io/).
 
-- User guide and tutorials: see the project documentation site
-- API reference: see the API section in the documentation site
+- User guide and tutorials: <https://graxpy.readthedocs.io/>
+- API reference: <https://graxpy.readthedocs.io/en/latest/api/index.html>
 
 For local docs builds from this repository, use:
 
@@ -31,14 +40,48 @@ names such as `"Si"` and `"Pt"`, with optional density overrides through
 and xrt-compatible material objects remain supported for now but emit a
 deprecation warning.
 
-PyPI project page: <https://pypi.org/project/graxpy/0.1.0/>
+PyPI project page: <https://pypi.org/project/graxpy/>
 
 For local editable installs:
 
 ```bash
 python -m pip install -e .
 ```
 
+## Select a solver
+
+RCWA is the default, so existing code continues to use it. Pass `solver` to
+choose the Nevière differential method instead:
+
+```python
+import grax
+
+grating = grax.LaminarGrating(
+    period_lpermm=400,
+    width_to_period_ratio=0.67,
+    depth_nm=14.9,
+    substrate_material="Si",
+    layer_material="Pt",
+    layer_thickness_nm=28.77,
+)
+
+common = dict(
+    grating=grating,
+    energy_ev=300.0,
+    grazing_angle_deg=4.0,
+    fourier_orders=30,
+    polarization="p",
+)
+
+rcwa_result = grax.run_simulation(**common, solver="rcwa")
+neviere_result = grax.run_simulation(**common, solver="neviere")
+```
+
+See the documentation's [Choosing a solver](docs/tutorials/choosing-a-solver.md)
+tutorial and the runnable
+`examples/simulation/neviere_solver/neviere_solver.py` for guidance and
+integration options.
+
 ## Local web app
 
 Install the package and web extra, then start the local server.
@@ -107,9 +150,9 @@ overlay for each run. Use `Manage runs` to rename runs or bulk delete them.
 
 ## Attribution
 
-`grax` is inspired by RETICOLO v9. This project is an independent Python
-implementation and is not an official RETICOLO distribution. RETICOLO is not
-bundled as part of the public `graxpy` package distribution.
+The modal RCWA solver is inspired by RETICOLO v9. GraxPy is an independent
+implementation, not an official RETICOLO port or distribution; RETICOLO is not
+bundled with the public `graxpy` package.
 
 - RETICOLO DOI: <https://doi.org/10.5281/zenodo.14631950>
 - RETICOLO license (CC BY 4.0): <https://creativecommons.org/licenses/by/4.0/>

diff --git a/docs/api/simulation/package-organization.md b/docs/api/simulation/package-organization.md
@@ -5,7 +5,7 @@ The `grax.simulation` module is implemented as a package to separate concerns by
 ## Module structure
 
 - `models.py`: Typed result dataclasses (`SingleSimulationResult`, `BatchSimulationResult`, `CaseExecutionResult`, etc.) and small compatibility containers.
-- `core.py`: One-point RCWA execution (`run_simulation`, `RCWASimulation`) and legacy plotting/comparison helpers.
+- `core.py`: One-point solver execution (`run_simulation`, `GratingSimulation`) and legacy plotting/comparison helpers.
 - `batch.py`: Generic batch execution (`BatchSimulationRunner`), checkpointing, subprocess execution, and worker calibration helpers.
 - `cases.py`: Lazy case-generation helpers (`energy_angle_cases`, `fixed_angle_cases`, `monochromator_cases`) and monochromator angle generation.
 - `theta_search.py`: Single-energy multilayer theta-search logic (`run_multilayer_theta_search`) and scan helpers.

diff --git a/docs/api/simulation/single-simulation.md b/docs/api/simulation/single-simulation.md
@@ -1,6 +1,7 @@
 # Single simulation
 
-Core API for one RCWA solve.
+Core API for one electromagnetic grating solve. Select the modal RCWA or
+Nevière differential-method path with `solver=`.
 
 ```{eval-rst}
 .. autofunction:: grax.run_simulation

diff --git a/docs/conf.py b/docs/conf.py
@@ -65,7 +65,7 @@
 }
 
 html_theme = "furo"
-html_static_path = ["_static"]
+html_static_path = ["_static"] if (Path(__file__).resolve().parent / "_static").is_dir() else []
 html_title = "grax documentation"
 templates_path = ["_templates"]
 

diff --git a/docs/developer/index.md b/docs/developer/index.md
@@ -1,12 +1,19 @@
 # Developer documentation
 
 This section is for contributors and maintainers who need to understand how the
-high-level API reaches the numerical RCWA implementation.
+high-level API reaches the numerical solvers.
+
+`grax` ships two selectable solver paths behind one interface: modal RCWA and
+the Nevière differential method. They share Fourier/discretization and
+boundary-response infrastructure, but use different layer-propagation methods.
+Read [RCWA theory](rcwa-theory.md) and the
+[Nevière differential method](neviere-theory.md) together for the two paths.
 
 ```{toctree}
 :maxdepth: 1
 
 rcwa-theory
+neviere-theory
 physics-to-code
 solver-flow
 module-guide

diff --git a/docs/developer/module-guide.md b/docs/developer/module-guide.md
@@ -14,17 +14,24 @@ This guide summarizes the source layout for contributors.
 - `stacks.py`: coating stack models above the substrate surface
 - `simulation/`: simulation workflows and orchestration (see package organization)
   - `models.py`: typed result dataclasses and small compatibility containers
-  - `core.py`: one-point RCWA execution and legacy plotting/comparison helpers
+  - `core.py`: one-point solver dispatch and legacy plotting/comparison helpers
   - `batch.py`: generic batch execution, checkpointing, subprocess execution, and worker calibration helpers
   - `cases.py`: lazy case-generation helpers and monochromator angle generation
   - `theta_search.py`: single-energy multilayer theta-search logic and scan helpers
   - `theta_search_sweep.py`: multi-energy adaptive theta-search sweep workflow and artifact writing
   - `serialization.py`: JSON/checkpoint serialization for single and case results
 - `parameter_sweep.py`: repeated simulations for numerical convergence and
   parameter studies
-- `rcwa_1d.py`: native one-dimensional TE-style RCWA numerical core
+- `solvers/`: the one-dimensional electromagnetic solvers
+  - `common.py`: shared types, `res0`/`res1`, the Fourier machinery, the layer
+    field operators (including the Li/fast-Fourier-factorization rules for TM),
+    the interface-response cascade, and the efficiency extraction
+  - `rcwa.py`: modal (eigen-decomposed) layer solve and `res2`
+  - `neviere.py`: Nevière differential method, `res2_dm`, and `NeviereOptions`
+- `rcwa_1d.py`: re-export shim keeping the historical `grax.rcwa_1d` import path
+  working. New code should import from `grax.solvers`
 - `roughness.py`: roughness post-processing models and diagnostics (currently
-  scalar Debye-Waller damping used by `rcwa_1d.res2()`)
+  scalar Debye-Waller damping applied by both solvers)
 - `slag.py`: legacy SLAG convenience helpers built on the generic public API
 
 ## Optimization package
@@ -48,15 +55,22 @@ to a scalar complex refractive index at a requested photon energy.
 Change `simulation/` when user-facing orchestration, result packaging,
 checkpointing, plotting, or validation behavior changes.
 
-Change `rcwa_1d.py` when the Fourier conversion, stack solve, polarization
-support, or diffraction-output calculation changes. Change `roughness.py` when
-roughness model math or diagnostics change.
+Change `solvers/common.py` when the Fourier conversion, the layer field
+operators, the interface cascade, or the diffraction-output calculation changes:
+that code is shared, so both solvers move together. Change `solvers/rcwa.py` or
+`solvers/neviere.py` only when one solver's own propagation changes. Adding a
+third solver means supplying a per-layer block builder to
+`solve_stack_from_layer_blocks()` and extending `SOLVER_NAMES` in
+`simulation/core.py`.
+
+Change `roughness.py` when roughness model math or diagnostics change.
 
 ## Reading order for new maintainers
 
-1. `examples/manual_single_simulation_example/manual_single_simulation.py`
+1. `examples/simulation/single_simulation/single_simulation.py`
 2. `examples/simulation/fixed_angle_sweep/fixed_angle_sweep.py`
 3. `src/grax/simulation/` (package overview in docs/api/simulation/package-organization.md)
 4. `src/grax/gratings.py`
 5. `src/grax/materials.py` and `src/grax/stacks.py`
-6. `src/grax/rcwa_1d.py`
+6. `src/grax/solvers/` (start with `common.py`, then read `rcwa.py` and
+   `neviere.py` as the two propagation implementations)