From 648744dfce7bd38ef81b4dc94da1561560915dc0 Mon Sep 17 00:00:00 2001 From: Robb Walters Date: Sat, 10 Oct 2026 11:01:45 -0700 Subject: [PATCH] feat(inelastic): refine band-bottom tables in lindhard run (#339) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The rate refinement now tests each cell at its quartiles and midpoint and, where those pass, at the vertex of the parabola through the largest sampled error and its two neighbours (parabolic interpolation of the peak). The midpoint alone (#351) missed in four tables because the error peaks off the midpoint; the quartiles alone still left the Al single pole at the se_yield.py grid at 1.0008 %. No constant is tuned: tolerance and window are #291's, the quartiles are the next bisection level, and the vertex is parameter-free. A cell exhausted to rounding is still an error. With it every refined table of the `rate` mode passes the dense 1 % gate (Al and Si, single pole and Mermin, 5 and 10 eV, 10 to 80 per decade), so `lindhard run` builds the inelastic table of every material with a band with RateRefinement::default(). The elastic table and the tables of bandless materials are unchanged. The table-cache key gains inelastic.rate_refinement (KEY_VERSION 4), a cached inelastic table must contain the run's grid in order, and electron_tables.csv matches elastic rows by energy. Committed δ(E) tables are stale pending #287. Closes #339 Loom-Issue: 2AMLogic/lindhard#339 Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 26 ++- docs/cli.md | 31 ++- docs/validation.md | 67 +++--- lindhard-cli/examples/inelastic_low_energy.rs | 3 +- lindhard-cli/src/electron.rs | 35 ++- lindhard-cli/src/table_cache.rs | 78 ++++++- lindhard-cli/tests/examples.rs | 44 ++++ lindhard/src/electron/inelastic/table.rs | 221 +++++++++++++++--- 8 files changed, 418 insertions(+), 87 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0aacb7e8..938fbfd5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,21 @@ one version). - `electron_psf_parameters.csv` gains `converged` (1/0) and `iterations` rows per fit, so the standalone CSV shows whether a fit hit its iteration budget (#347). +### Changed + +- `lindhard run` builds the inelastic table of every material with a `band` + with the rate refinement (`RateRefinement::default()`, #339): rows are + added between 2 and 20 eV above the Fermi level until the interpolated + inelastic rate is within 1 % of the exact one. **Results of runs with a + `band` change**; the elastic table and the tables of materials without a + band are unchanged, bit for bit. The committed secondary-yield δ(E) tables + (`validation/experiments/se_yield_results.json`) predate it and are stale + until their re-run (#287). `electron_tables.csv` has a line per inelastic + row, with an empty elastic column at the added rows (the elastic value is + now matched by energy). The table-cache key gains `inelastic.rate_refinement` + (`KEY_VERSION` 4), so cached tables are rebuilt once, and a cached + inelastic table must contain the run's grid in order rather than equal it. + ### Fixed - Python `RunResult.write()` now follows the command's output lifecycle: it @@ -32,16 +47,17 @@ one version). of it): rows added to one band-bottom inelastic table from the Fermi level upward. Every grid cell overlapping `E - E_F` = 2 to 20 eV is bisected at its midpoint while the linear interpolation of the exact rate - misses it there by more than 1 %; more than `max_rows` added rows is an + misses it by more than 1 % at the cell's quartiles, midpoint or estimated + error peak (the vertex of the parabola through the largest sampled error + and its neighbours); more than `max_rows` added rows is an error, and so is a failing cell bisected down to rounding (a step in the - rate), each naming the cell. The grid is the same on any thread count. Off by default, so every - existing table is unchanged, and refused on the model's own axis. + rate), each naming the cell. The grid is the same on any thread count. Off by default, and + refused on the model's own axis. `RATE_TOLERANCE` and `RATE_WINDOW_EV` (the #291 gate, moved from the example into `electron::inelastic::table`) are its defaults. The `rate` mode of `lindhard-cli/examples/inelastic_low_energy.rs` also evaluates each coarse table refined (`docs/validation.md`, "Table grid above the - Fermi level"). `lindhard run` does not use it yet: the Al single pole at - the `se_yield.py` grid misses the 1 % gate by 0.0008 points. + Fermi level"); every refined table there passes the 1 % gate. - `[electron.materials..mermin_fit]` (#306): per-material options of the `mermin-melf` oscillator fit, `oscillators` (1 to 16, default 3), `weighting` (`"relative"`, the default, or `"uniform"`) and diff --git a/docs/cli.md b/docs/cli.md index b30736a5..341d01bf 100644 --- a/docs/cli.md +++ b/docs/cli.md @@ -512,7 +512,9 @@ fails, naming the key: remove the key. Results of runs with a `band` change (see `CHANGELOG.md`). **`[electron.tables]`**: one log-spaced energy grid shared by the elastic and -inelastic tables of every material. +inelastic tables of every material. The inelastic table of a material with +a `band` holds extra rows besides (#339; "Grid for a material with a band", +below). | Key | Default | Meaning | |---|---|---| @@ -539,11 +541,15 @@ grid above the Fermi level"; reproduce with can add rows to the inelastic table from `E_F` upward (`InelasticTableOptions::rate_refinement`, #339), which brings every one of those tables to about 1 % at both the default grid and the grid of -`validation/experiments/se_yield.py`, but one of them (the Al single pole -from 5 eV at 10 points per decade) misses 1 % by 0.0008 points, so -`lindhard run` does not use it yet (#339). Until then a denser grid helps -but does not close the gap: 40 points per decade cuts the error by a -factor of 2 to 4 from 20. +`validation/experiments/se_yield.py`. **`lindhard run` uses it for the +inelastic table of every material with a `band`** (#339): every cell of the +grid overlapping `E - E_F` = 2 to 20 eV is bisected while the linear +interpolation of the rate misses the exact rate by more than 1 % at the +cell's quartiles, midpoint or estimated error peak. The table then has the +grid's rows and the added ones (a few tens: Al and Si, 5 to 25 rows at +these grids); `physics.materials.inelastic_table.energies` counts them. The +elastic table, and the inelastic table of a material without a band, keep +exactly the grid's rows. A denser grid is not needed for this region. **`[electron.materials.]`**: the electron data of each material the target uses, keyed by the name the target gives it (a `[materials]` key or an @@ -816,8 +822,10 @@ skip them. Parameters from a fit with `converged` 0 need inspection of the fit diagnostics before use. `electron_tables.csv`: `material,energy_ev,elastic_inverse_mfp_per_nm,inelastic_inverse_mfp_per_nm,inelastic_mean_loss_ev,inelastic_stopping_ev_per_nm`, -the tables the run used, per material and grid energy (the stopping power is -`λ⁻¹ ⟨W⟩` of the stored loss distribution). +the tables the run used, per material and inelastic-table energy (the +stopping power is `λ⁻¹ ⟨W⟩` of the stored loss distribution). The elastic +column is empty at the rows the rate refinement adds to the inelastic table +of a material with a band (#339), which the elastic table does not have. ### Cross-section table cache (`--table-cache`) @@ -851,7 +859,9 @@ the table depends on: SHA-256 and provenance of the optical ELF file, the material's band parameters (they set the table's axis and Fermi energy for a material with a band) and, with `mermin-melf`, the material's resolved Mermin fit options - (`mermin_fit`; no table and a table of the defaults give the same key). + (`mermin_fit`; no table and a table of the defaults give the same key), + and the rate refinement of a material with a band (`rate_refinement`: + its window, tolerance and row bound; `null` without a band). Every `f64` is written in shortest round-trip form, so a change in the last bit of any number is a different key. A lookup must find the stored key @@ -862,7 +872,8 @@ corrupted anywhere, even in a single cross-section or probability value that still parses, is refused with an error naming the file, before it is parsed, and a missing `.sha256` file is an error too. The table is then loaded and validated by the library, and its axis and energy grid are checked against the -run. A file found under the run's key that fails any of these checks is an +run: an elastic table must have exactly the run's grid, an inelastic table +must contain it in order (it may hold the refinement's rows as well). A file found under the run's key that fails any of these checks is an error, never a silent rebuild; remove the entry's files to rebuild it. A table of another cache `format_version` can never be found, since the version is part of the key. Writes go to a temporary file renamed into place, the table diff --git a/docs/validation.md b/docs/validation.md index d6eb2aa5..76b6226e 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -1994,38 +1994,51 @@ result worse, not better; no Cu, Si or Au set changes status. rate needs rows from `E_F` upward, not a denser uniform grid. **Rate refinement (#339).** `InelasticTableOptions::rate_refinement` - (`electron::inelastic::table`, "Rate refinement"; off by default) gives - one band-bottom table those rows: every grid cell overlapping `E - E_F` - = 2 to 20 eV is bisected at its midpoint while the linear interpolation - of the exact rate misses the exact rate there by more than 1 % - (`RATE_TOLERANCE` and `RATE_WINDOW_EV`, the gate's own constants, now in - the library and read by the example). The `rate` mode builds each coarse - table again with it (`+refined`: the coarse rows and every row the - construction adds) and gates it the same way. Largest error over 2 to 20 - eV (where it falls) and rows of the refined table (the coarse rows + (`electron::inelastic::table`, "Rate refinement") gives one band-bottom + table those rows: every grid cell overlapping `E - E_F` = 2 to 20 eV is + bisected at its midpoint while the linear interpolation of the exact + rate misses the exact rate by more than 1 % (`RATE_TOLERANCE` and + `RATE_WINDOW_EV`, the gate's own constants, in the library and read by + the example) at any of the cell's test points: its quartiles and + midpoint, and, where those pass, the vertex of the parabola through the + largest sampled error and its two neighbours (the estimated peak). The + `rate` mode builds each coarse table again with it (`+refined`: the + coarse rows and every row the construction adds) and gates it the same + way. The first version (#351) tested the midpoint only and missed in + four tables, by 1.0008 % (Al single pole, 5 eV, 10 per decade, at 6.1 + eV) to 1.43 %: there the error peaks off the midpoint. Testing the + quartiles as well still left the first of them at 1.0008 % (its peak lies + between the first quartile and the midpoint of the cell 5.3 to 7.2 eV, + which miss by 0.84 and 0.99 %) and Al Mermin from 10 eV at 10 per decade + at 1.03 %; the estimated peak closes both. No constant was tuned: the + tolerance and window are #291's, and the test points are the next level + of the bisection and the parabolic vertex. Largest error over 2 to 20 eV + (where it falls) and rows of the refined table (the coarse rows bracketing 0.5 to 20 eV above `E_F` plus the rows added), at the two grids that matter, the `se_yield.py` grid (5 eV, 10 per decade) and the - library default (10 eV, 20 per decade): + library default (10 eV, 20 per decade), 2026-10-10: | Table | 5 eV, 10 / decade | rows (coarse + added) | 10 eV, 20 / decade | rows (coarse + added) | |---|---|---|---|---| - | Al single-pole | **1.0008 % (6.1), fails** | 6 + 22 | 0.84 % (3.1) | 11 + 24 | - | Al Mermin | 0.66 % (18.6) | 6 + 10 | 1.00 % (2.7) | 11 + 6 | - | Si single-pole | 0.84 % (15.2) | 6 + 25 | 0.96 % (3.9) | 10 + 21 | - | Si Mermin | 0.93 % (3.0) | 6 + 7 | 0.72 % (15.8) | 10 + 5 | - - At the other grids of the mode (5 or 10 eV, 10 to 80 per decade) the - refined tables pass except the Al single pole from 5 eV at 20 per - decade (1.43 % at 17.6 eV), Al Mermin from 10 eV at 10 per decade - (1.03 %, 4.6 eV) and Si Mermin from 10 eV at 10 per decade (1.40 %, - 18.6 eV); the uniform grids alone fail everywhere below 80 per decade. - The midpoint test is a criterion, not a bound: in those cells the error - peaks away from the midpoint the test samples. By the stop rule of #339 - the criterion is not tightened by an empirical factor and neither the - tolerance nor the window is widened, so **`lindhard run` does not use - the refinement yet** and the grid of these runs is unchanged. #339 stays - open for that decision; #291 and #287 are not settled by it. - The committed results are unchanged. + | Al single-pole | 0.967 % (16.4) | 6 + 23 | 0.843 % (3.1) | 11 + 24 | + | Al Mermin | 0.656 % (18.6) | 6 + 10 | 0.880 % (4.8) | 11 + 7 | + | Si single-pole | 0.843 % (15.2) | 6 + 25 | 0.953 % (4.9) | 10 + 22 | + | Si Mermin | 0.934 % (3.0) | 6 + 7 | 0.721 % (15.8) | 10 + 5 | + + Every refined table of the mode passes (5 and 10 eV, 10 to 80 per + decade, 32 tables), including the three that failed with the midpoint + test outside these grids: the Al single pole from 5 eV at 20 per decade + 0.772 % (15.9; 11 + 23 rows), Al Mermin from 10 eV at 10 per decade + 0.893 % (2.7; 7 + 11) and Si Mermin from 10 eV at 10 per decade 0.842 % + (18.2; 6 + 10). The uniform grids alone fail everywhere below 80 per + decade. The test is a criterion, not a bound; the dense gate is the + check. **`lindhard run` now builds the inelastic table of every material + with a band with the refinement** (`RateRefinement::default()`; the + elastic table and tables without a band keep the grid), which settles the + grid question of #291 for the δ runs: the `se_yield.py` grid (5 eV, 10 + per decade) needs no change. The committed δ tables below were computed + before it, so they are stale until the re-run of #287; this change does + not re-run them. - Secondary generation: Kieft and Bosch (2008) as implemented (`secondaries = "kieft-bosch"`); transport cutoff 1 eV above the vacuum level (`cutoff_reference = "vacuum-level"`) unless a row says otherwise. diff --git a/lindhard-cli/examples/inelastic_low_energy.rs b/lindhard-cli/examples/inelastic_low_energy.rs index 1dbaf3f0..550b0b54 100644 --- a/lindhard-cli/examples/inelastic_low_energy.rs +++ b/lindhard-cli/examples/inelastic_low_energy.rs @@ -53,7 +53,8 @@ //! builder refuses is printed as refused. Each coarse table is followed by //! the same rows built with the library's rate refinement //! (`InelasticTableOptions::rate_refinement`, `RateRefinement::default()`, -//! #339; `+refined`; `lindhard run` does not use it yet): its grid is the +//! #339; `+refined`; `lindhard run` uses it for every table on the +//! band-bottom axis): its grid is the //! coarse rows and every row the refinement adds, and `rows` counts them. //! The largest error is printed to three decimals, so that a verdict at //! the tolerance can be read. The verdict that 20 and 40 points per decade fail diff --git a/lindhard-cli/src/electron.rs b/lindhard-cli/src/electron.rs index dbd5d4cb..614d492b 100644 --- a/lindhard-cli/src/electron.rs +++ b/lindhard-cli/src/electron.rs @@ -22,6 +22,7 @@ use lindhard::electron::elastic::table::{ use lindhard::electron::elastic::{SalvatDhfs, SolverOptions}; use lindhard::electron::inelastic::table::{ build_inelastic_table_for_model, mean_loss_ev, EnergyAxis, InelasticTableOptions, + RateRefinement, }; use lindhard::electron::inelastic::{PennAlgorithm, PennInelastic}; use lindhard::electron::transport::{LayerTables, RunMetadata, Transport}; @@ -236,13 +237,30 @@ fn inelastic_table( other => PennInelastic::try_new(other, m.optical_elf.clone())?, } .with_fermi_energy_ev(fermi_ev)?; + let mut options = InelasticTableOptions::new(r.table_energy_ev.clone()).with_axis(axis); + options.rate_refinement = rate_refinement(r, m); Ok(build_inelastic_table_for_model( &model, &m.material, - &InelasticTableOptions::new(r.table_energy_ev.clone()).with_axis(axis), + &options, )?) } +/// The rate refinement of the inelastic table of `m` (#339): the library's +/// default ([`RateRefinement::default`]) on the band-bottom axis of a +/// material with a band, where the rate climbs from zero at the Fermi +/// level; `None` on the model's own axis, whose table keeps exactly the +/// rows of `electron.tables`. The refined table has the grid's rows and the +/// rows the refinement adds between 2 and 20 eV above the Fermi level +/// (`lindhard::electron::inelastic::table`, "Rate refinement"); the elastic +/// table keeps the grid. Part of the table cache key. +pub fn rate_refinement( + r: &ResolvedElectron, + m: &ResolvedElectronMaterial, +) -> Option { + (inelastic_axis(r, m).1 == EnergyAxis::BandBottom).then(RateRefinement::default) +} + /// The model's Fermi energy and the table axis of the inelastic table of `m`. /// /// For a material with a band, the table is built on the band-bottom axis @@ -740,14 +758,15 @@ pub fn tables_csv(r: &ResolvedElectron, sim: &ElectronSimulation) -> String { inelastic_mean_loss_ev,inelastic_stopping_ev_per_nm\n", ); for (m, t) in r.materials.iter().zip(&sim.tables) { - // Both tables share the grid of electron.tables. + // One line per inelastic row. The elastic table has the rows of + // electron.tables, which the inelastic table contains; the rows the + // rate refinement adds (#339) have no elastic value. + let elastic = t.elastic.energy_ev(); for (i, &e) in t.inelastic.energy_ev().iter().enumerate() { - let el = t - .elastic - .energy_ev() - .get(i) - .filter(|&&x| x == e) - .map(|_| t.elastic.inverse_mfp_per_m()[i] * NM); + let el = elastic + .binary_search_by(|x| x.total_cmp(&e)) + .ok() + .map(|j| t.elastic.inverse_mfp_per_m()[j] * NM); let inv = t.inelastic.inverse_mfp_per_m()[i]; let mean = t .inelastic diff --git a/lindhard-cli/src/table_cache.rs b/lindhard-cli/src/table_cache.rs index 9b288209..161d1f23 100644 --- a/lindhard-cli/src/table_cache.rs +++ b/lindhard-cli/src/table_cache.rs @@ -36,14 +36,17 @@ //! material (the table of a material with a band is built on the //! band-bottom axis with the band's Fermi energy, `electron::inelastic_axis`), //! and, for `mermin-melf`, the options of the material's oscillator fit -//! (`[electron.materials..mermin_fit]`, #306). +//! (`[electron.materials..mermin_fit]`, #306), and the rate +//! refinement of a band-bottom table (`electron::rate_refinement`, #339). //! //! A cache file of another format version can never be found under a //! current key (the version is in the key). A file that **is** found under //! the run's key is loaded with the library loader, whose errors (an //! unsupported version, a failed validation) are returned, not treated as a //! miss; after loading, its sampling axis and energy grid are checked -//! against the run. +//! against the run: an elastic table must have exactly the run's grid, an +//! inelastic table must contain it in order (it may hold the rows of its +//! rate refinement as well). //! //! Writes go to a temporary file in the same directory and are renamed into //! place, the table first, then its hash, then its key, so the key is the @@ -65,7 +68,9 @@ use sha2::{Digest, Sha256}; /// removed or changes meaning. /// /// 3: `inelastic.mermin_fit`, the Mermin fit options (#306). -pub const KEY_VERSION: u32 = 3; +/// 4: `inelastic.rate_refinement`, the rate refinement of a band-bottom +/// table (#339). +pub const KEY_VERSION: u32 = 4; /// Lowercase hex SHA-256 of `bytes`. pub fn sha256_hex(bytes: &[u8]) -> String { @@ -159,6 +164,11 @@ struct InelasticKey { /// [`lindhard::electron::inelastic::MerminFitOptions`] with /// `mermin-melf`; `null` with the other models, which do not fit. mermin_fit: Option, + /// `Debug` of the table's rate refinement + /// ([`crate::electron::rate_refinement`]; `null` without one): its + /// window, tolerance and row bound. The rows it adds are a deterministic + /// function of these, the model and the grid. + rate_refinement: Option, } #[derive(Serialize)] @@ -208,6 +218,7 @@ pub fn key_text( band: format!("{:?}", m.band), mermin_fit: (r.inelastic == PennAlgorithm::Mermin) .then(|| format!("{:?}", m.mermin_fit)), + rate_refinement: crate::electron::rate_refinement(r, m).map(|x| format!("{x:?}")), }), ), }; @@ -461,11 +472,22 @@ impl TableCache { kind.axis() ); } - if table.energy_ev() != r.table_energy_ev.as_slice() { + // The elastic table has exactly the run's grid. The inelastic table + // has the run's grid and, on the band-bottom axis, the rows of its + // rate refinement (#339), so it must contain the grid in order. + let grid_ok = match kind { + TableKind::Elastic => table.energy_ev() == r.table_energy_ev.as_slice(), + TableKind::Inelastic => contains_in_order(table.energy_ev(), &r.table_energy_ev), + }; + if !grid_ok { bail!( - "table cache {}: field `energy_ev` differs from this run's grid \ + "table cache {}: field `energy_ev` {} this run's grid \ ({} energies cached, {} needed)", table_path.display(), + match kind { + TableKind::Elastic => "differs from", + TableKind::Inelastic => "does not contain", + }, table.energy_ev().len(), r.table_energy_ev.len() ); @@ -482,6 +504,13 @@ impl TableCache { } } +/// Whether `grid` is a subsequence of `cached` (every energy of `grid`, in +/// order, bit for bit). +fn contains_in_order(cached: &[f64], grid: &[f64]) -> bool { + let mut it = cached.iter(); + grid.iter().all(|g| it.any(|c| c.to_bits() == g.to_bits())) +} + #[cfg(test)] mod tests { use super::*; @@ -698,6 +727,33 @@ mod tests { }); } + /// A band-bottom table is refined (#339) and keyed so; a table without + /// a band is not, and its key says `null`. + #[test] + fn key_carries_the_rate_refinement() { + let r = example(); + let (_, i) = keys(&r); + assert!(i.contains("\"rate_refinement\":\"RateRefinement {"), "{i}"); + let mut bandless = r.clone(); + bandless.materials[0].band = None; + assert!(crate::electron::rate_refinement(&bandless, &bandless.materials[0]).is_none()); + let (_, i) = keys(&bandless); + assert!(i.contains("\"rate_refinement\":null"), "{i}"); + } + + #[test] + fn grid_containment_is_in_order_and_bit_exact() { + let grid = [1.0, 2.0, 4.0]; + assert!(contains_in_order(&grid, &grid)); + assert!(contains_in_order(&[1.0, 1.5, 2.0, 3.0, 4.0, 8.0], &grid)); + assert!(!contains_in_order(&[1.0, 4.0, 2.0], &grid)); + assert!(!contains_in_order(&[1.0, 2.0], &grid)); + assert!(!contains_in_order( + &[1.0, f64::from_bits(2.0f64.to_bits() + 1), 4.0], + &grid + )); + } + #[test] fn key_changes_with_the_build() { let r = example(); @@ -820,6 +876,18 @@ mod tests { "\"band\":{}", serde_json::to_string(&format!("{:?}", r.materials[0].band)).unwrap() ))); + // The band also sets the rate refinement of the inelastic table + // (`electron::rate_refinement`, #339), the function the build uses: + // `inelastic.rate_refinement`. + let refinement = crate::electron::rate_refinement(&r, &r.materials[0]); + assert!(refinement.is_some()); + assert!( + i.contains(&format!( + "\"rate_refinement\":{}", + serde_json::to_string(&format!("{:?}", refinement.unwrap())).unwrap() + )), + "{i}" + ); let MerminFitSpec { // `inelastic.mermin_fit` (`n_oscillators`). oscillators, diff --git a/lindhard-cli/tests/examples.rs b/lindhard-cli/tests/examples.rs index 5085aeee..e218bc2a 100644 --- a/lindhard-cli/tests/examples.rs +++ b/lindhard-cli/tests/examples.rs @@ -1402,6 +1402,20 @@ fn electron_table_cache_reuse_is_bit_identical() { ); let plain = json(&out("plain").join("electron_summary.json")); assert!(plain["physics"]["materials"][0]["inelastic_table"]["cache"].is_null()); + // The example's material has a band, so its inelastic table holds the + // rows of the rate refinement (#339) as well: a refined table read back + // from the cache gives the same results, bit for bit (below). + let energies = |t: &str| { + plain["physics"]["materials"][0][t]["energies"] + .as_u64() + .unwrap() + }; + assert!( + energies("inelastic_table") > energies("elastic_table"), + "{} vs {}", + energies("inelastic_table"), + energies("elastic_table") + ); let miss = json(&out("miss").join("electron_summary.json")); let hit = json(&out("hit-1").join("electron_summary.json")); for t in ["elastic_table", "inelastic_table"] { @@ -1486,6 +1500,36 @@ fn electron_table_cache_misses_on_changed_physics_and_refuses_bad_files() { ok(&go(&src, "a")); assert_eq!(cache_files(&cache).len(), 6); + // The inelastic table of the example's band material holds the rows of + // its rate refinement (#339) besides the run's grid. One whose grid no + // longer contains the run's grid (an energy of the grid edited, the + // hash with it) is refused; restored, it is read. + let ikey = cache_files(&cache) + .into_iter() + .find(|f| f.starts_with("inelastic-") && f.ends_with(".key.json")) + .unwrap(); + let itable_path = cache.join(ikey.replace(".key.json", ".toml")); + let itable = std::fs::read_to_string(&itable_path).unwrap(); + let first = itable.find("energy_ev = [").unwrap() + "energy_ev = [".len(); + let end = first + itable[first..].find(',').unwrap(); + write_table_and_hash( + &itable_path, + &format!("{}10.5{}", &itable[..first], &itable[end..]), + ); + let o = go(&src, "i"); + assert!(!o.status.success()); + let e = String::from_utf8_lossy(&o.stderr); + assert!( + e.contains("field `energy_ev` does not contain this run's grid"), + "{e}" + ); + write_table_and_hash(&itable_path, &itable); + ok(&go(&src, "j")); + assert_eq!( + table_sources(&dir.join("j")), + ("cache".into(), "cache".into()) + ); + // Another band with the same inner potential (so the same table grid): // the inelastic table, built on the band-bottom axis with the band's // minimum excitation energy (#241), is built and stored beside the diff --git a/lindhard/src/electron/inelastic/table.rs b/lindhard/src/electron/inelastic/table.rs index dabd66cc..114e8945 100644 --- a/lindhard/src/electron/inelastic/table.rs +++ b/lindhard/src/electron/inelastic/table.rs @@ -158,13 +158,17 @@ //! every cell of the grid that overlaps `E - E_F` in //! [`RateRefinement::window_ev`] ([`RATE_WINDOW_EV`] by default) is bisected //! at its midpoint in `E` while the linear interpolation of the exact row -//! rate between the cell's ends differs from the exact rate at the midpoint -//! by more than [`RateRefinement::tolerance`] ([`RATE_TOLERANCE`]) of it. +//! rate between the cell's ends differs from the exact rate at any of the +//! cell's test points by more than [`RateRefinement::tolerance`] +//! ([`RATE_TOLERANCE`]) of it. The test points are the cell's quartiles and +//! midpoint (`rate_test_points`) and, where those pass, the estimated +//! position of the largest error: the vertex of the parabola through the +//! largest of the sampled errors and its two neighbours (`rate_error_peak`). //! Children that no longer overlap the window are not tested, so the //! bisection never approaches the zero of the rate at `E_F`, and cells //! outside the window are never split. Only the scalar rate is evaluated //! while refining (the row rule of the table: zero for an empty row); each -//! level's midpoints are evaluated in parallel and merged in sorted order, +//! level's test points are evaluated in parallel and merged in sorted order, //! so the grid is the same on any thread count. The rows are then built as //! any other, from the merged grid, which contains the given grid. A //! refinement that would add more than [`RateRefinement::max_rows`] rows is @@ -175,9 +179,13 @@ //! rounding (about 50 levels) before any useful row bound: that second //! error, not the row bound, is what stops it. //! -//! The midpoint test is a criterion, not a bound: for a smooth rate the -//! midpoint error of a cell is close to its largest, but the acceptance -//! check is the dense evaluation of the `rate` mode of +//! The test is a criterion, not a bound. A test of the midpoint alone +//! (#351) left cells whose error peaked off the midpoint above the +//! tolerance, and so did one of the quartiles and midpoint (the Al single +//! pole from 5 eV at 10 per decade, 1.0008 % at 6.1 eV, between the first +//! quartile and the midpoint of its cell); the parabolic vertex locates such +//! a peak without a tuned number of samples. The acceptance check is the dense +//! evaluation of the `rate` mode of //! `lindhard-cli/examples/inelastic_low_energy.rs` (every 0.1 eV over the //! window, `docs/validation.md`). The table's other grid (the elastic table //! of a run, or a table without the option) is unchanged: the transport @@ -569,7 +577,7 @@ pub const RATE_TOLERANCE: f64 = 0.01; pub const RATE_WINDOW_EV: (f64, f64) = (2.0, 20.0); /// The default bound on the rows the rate refinement may add (a safety -/// bound against a rate that fails the midpoint test in many cells at once; +/// bound against a rate that fails the test in many cells at once; /// not a tuned value): 1000 rows, so cells of 18 meV on average across the /// 18 eV of [`RATE_WINDOW_EV`], far below the 0.1 eV step of the `rate` /// gate. A step in the rate splits one cell per level and is stopped @@ -585,7 +593,7 @@ pub struct RateRefinement { /// refined when it overlaps it. pub window_ev: (f64, f64), /// The relative tolerance of the linear interpolation of the rate at a - /// cell's midpoint. + /// cell's test points (its quartiles and midpoint). pub tolerance: f64, /// The most rows the refinement may add; more is an error. pub max_rows: usize, @@ -972,6 +980,64 @@ fn row_rate(penn: &M, axis: EnergyAxis, t: f64) -> Result 0.0 { inv } else { 0.0 }) } +/// The points at which the rate refinement tests the cell `a..b` (module +/// docs, "Rate refinement"): its first quartile, midpoint and third +/// quartile, each formed as the midpoint of a half, so that a quartile is +/// bitwise the midpoint of the half a bisection would make. +/// +/// The interpolation error of a function with constant curvature peaks at +/// the midpoint; the band-bottom rate's curvature changes across a cell, +/// and the relative error divides by a rate that changes across it, so the +/// peak moves off the midpoint (#339 measured misses of up to 1.43 % in +/// cells whose midpoint missed by less than 1 %). The quartiles are the +/// points of the next level of the bisection itself, a fixed dyadic set +/// rather than a tuned count; they bracket the peak, which +/// `rate_error_peak` then locates. +fn rate_test_points(a: f64, b: f64) -> [f64; 3] { + let m = 0.5 * (a + b); + [0.5 * (a + m), m, 0.5 * (m + b)] +} + +/// The estimated position of the largest relative interpolation error in +/// the cell `a..b`, from its signed relative errors `e` at the quartiles and +/// midpoint (`rate_test_points`; the error is zero at the cell's ends): the +/// vertex of the parabola through the largest of the five samples and its +/// two neighbours (parabolic interpolation of an extremum, the step of +/// Brent's method; R. P. Brent, *Algorithms for Minimization without +/// Derivatives*, Prentice-Hall 1973, ch. 5). The largest sample's +/// neighbours do not exceed it in magnitude, so the vertex lies between +/// them. `None` where every sample is zero or the three are collinear. +fn rate_error_peak(a: f64, b: f64, e: [f64; 3]) -> Option { + let s = [0.0, e[0], e[1], e[2], 0.0]; + // The largest sample in magnitude (the midpoint on a tie). + let k = [2, 1, 3] + .into_iter() + .fold(2, |k, j| if s[j].abs() > s[k].abs() { j } else { k }); + let (l, c, r) = (s[k - 1], s[k], s[k + 1]); + let curvature = l - 2.0 * c + r; + if c == 0.0 || curvature == 0.0 || !curvature.is_finite() { + return None; + } + let h = 0.25 * (b - a); + let xk = a + h * k as f64; + let v = xk + 0.5 * h * (l - r) / curvature; + (v > a && v < b && v.is_finite()).then_some(v) +} + +/// The miss of the interpolated rate at a test point of the rate refinement. +#[derive(Debug, Clone, Copy)] +struct RateMiss { + /// `(interpolated - exact) / exact` (infinite where the exact rate is + /// zero and the interpolated one is not). + rel: f64, + /// The energy of the test point, eV. + at: f64, + /// The interpolated rate, 1/m. + interpolated: f64, + /// The exact rate, 1/m. + exact: f64, +} + /// The energy grid of a table with the rows of the rate refinement of /// `options` merged in (module docs, "Rate refinement"), or `None` without /// one. `options` is checked. @@ -1017,7 +1083,7 @@ fn refined_energy_grid( "rate refinement", format!( "cannot hold the rate to {} between E - E_F = {} and {} eV: cell {}..{} eV \ - is exhausted to rounding and its parent still failed at its midpoint (end \ + is exhausted to rounding and its parent still failed at a test point (end \ rates {:e} and {:e} 1/m)", r.tolerance, r.window_ev.0, @@ -1032,19 +1098,72 @@ fn refined_energy_grid( if cells.is_empty() { break; } - let mids: Vec = cells.iter().map(|&i| 0.5 * (x[i] + x[i + 1])).collect(); - let exact: Vec = mids + // The test points of each cell, in cell order: its quartiles and + // midpoint (`rate_test_points`), then the estimated peak of its + // error (`rate_error_peak`) where those pass. + let points: Vec<[f64; 3]> = cells + .iter() + .map(|&i| rate_test_points(x[i], x[i + 1])) + .collect(); + let exact: Vec = points .par_iter() - .map(|&e| rate(e)) + .flat_map_iter(|p| p.iter().copied()) + .map(rate) .collect::>()?; - // Split where the transport's linear interpolation misses the exact - // rate at the midpoint by more than the tolerance (a zero exact rate - // with a nonzero interpolant always fails). - let split: Vec = cells + // The miss of the transport's linear interpolation in cell `i` at + // `e`, where the exact rate is `m`: signed, relative to `m` (a zero + // exact rate with a nonzero interpolant is an infinite miss), with + // the interpolated rate. + let miss = |i: usize, e: f64, m: f64| { + let lin = y[i] + (e - x[i]) / (x[i + 1] - x[i]) * (y[i + 1] - y[i]); + let d = lin - m; + let rel = if d == 0.0 { + 0.0 + } else if m > 0.0 { + d / m + } else { + f64::INFINITY + }; + RateMiss { + rel, + at: e, + interpolated: lin, + exact: m, + } + }; + let mut worst: Vec = Vec::with_capacity(cells.len()); + let mut peaks: Vec<(usize, f64)> = Vec::new(); + for (k, (&i, (p, ex))) in cells .iter() - .zip(&exact) - .map(|(&i, &m)| (0.5 * (y[i] + y[i + 1]) - m).abs() > r.tolerance * m) - .collect(); + .zip(points.iter().zip(exact.chunks_exact(3))) + .enumerate() + { + let m: Vec = p.iter().zip(ex).map(|(&e, &v)| miss(i, e, v)).collect(); + let w = m + .iter() + .fold(m[1], |a, &b| if b.rel.abs() > a.rel.abs() { b } else { a }); + if w.rel.abs() <= r.tolerance { + if let Some(e) = rate_error_peak(x[i], x[i + 1], [m[0].rel, m[1].rel, m[2].rel]) { + peaks.push((k, e)); + } + } + worst.push(w); + } + // The estimated peaks of the cells that passed at their quartiles and + // midpoint, in parallel and in cell order. + let peak_exact: Vec = peaks + .par_iter() + .map(|&(_, e)| rate(e)) + .collect::>()?; + for (&(k, e), &v) in peaks.iter().zip(&peak_exact) { + let p = miss(cells[k], e, v); + if p.rel.abs() > worst[k].rel.abs() { + worst[k] = p; + } + } + // Split at the midpoint where the miss at any test point exceeds the + // tolerance. + let split: Vec = worst.iter().map(|w| w.rel.abs() > r.tolerance).collect(); let n_split = split.iter().filter(|&&s| s).count(); if n_split == 0 { break; @@ -1052,20 +1171,27 @@ fn refined_energy_grid( if added + n_split > r.max_rows { let k = split.iter().position(|&s| s).expect("one split"); let i = cells[k]; + let RateMiss { + at, + interpolated: lin, + exact: m, + .. + } = worst[k]; return Err(invalid( "rate refinement", format!( "more than {} rows needed to hold the rate to {} between E - E_F = {} and \ - {} eV (cell {}..{} eV still fails at its midpoint: interpolated {:e}, \ - exact {:e} 1/m)", + {} eV (cell {}..{} eV still fails at {} eV: interpolated {:e}, exact {:e} \ + 1/m)", r.max_rows, r.tolerance, r.window_ev.0, r.window_ev.1, x[i], x[i + 1], - 0.5 * (y[i] + y[i + 1]), - exact[k] + at, + lin, + m ), )); } @@ -1079,9 +1205,10 @@ fn refined_energy_grid( ny.push(y[i]); if k < cells.len() && cells[k] == i { if split[k] { - let m = mids[k]; + // The midpoint, a test point: its exact rate is known. + let m = points[k][1]; nx.push(m); - ny.push(exact[k]); + ny.push(exact[3 * k + 1]); ntodo.push(overlaps(x[i], m)); ntodo.push(overlaps(m, x[i + 1])); } else { @@ -1129,7 +1256,8 @@ fn refinement_provenance(options: &InelasticTableOptions, added: usize) -> Strin None => String::new(), Some(r) => format!( "; rate refinement (#339): {added} rows added where the linear interpolation of \ - the rate missed it by more than {} at a cell midpoint, E - E_F = {}..{} eV", + the rate missed it by more than {} at a cell's quartiles, midpoint or estimated \ + error peak, E - E_F = {}..{} eV", r.tolerance, r.window_ev.0, r.window_ev.1 ), } @@ -1786,7 +1914,7 @@ mod tests { } #[test] - fn rate_refinement_adds_rows_only_in_the_window_and_meets_the_midpoint_test() { + fn rate_refinement_adds_rows_only_in_the_window_and_meets_the_test() { let model = SteepRate(steep); let grid = vec![5.0, 9.0, 14.0, 22.0, 35.0, 60.0, 100.0]; let o = band_options(grid.clone()); @@ -1804,13 +1932,16 @@ mod tests { let (a, b) = (grid[k - 1], grid[k]); assert!(b > lo && a < hi, "{e} in the cell {a}..{b}"); } - // Every final cell overlapping the window passes the midpoint test. + // Every final cell overlapping the window passes the test at its + // quartiles and midpoint. let rate = |e: f64| row_rate(&model, EnergyAxis::BandBottom, e - 10.0).unwrap(); for w in refined.windows(2) { if w[1] > lo && w[0] < hi { - let m = rate(0.5 * (w[0] + w[1])); - let lin = 0.5 * (rate(w[0]) + rate(w[1])); - assert!((lin - m).abs() <= RATE_TOLERANCE * m, "{w:?}"); + for e in rate_test_points(w[0], w[1]) { + let m = rate(e); + let lin = rate(w[0]) + (e - w[0]) / (w[1] - w[0]) * (rate(w[1]) - rate(w[0])); + assert!((lin - m).abs() <= RATE_TOLERANCE * m, "{w:?} at {e}"); + } } } // Without the option the grid is the given one. @@ -1882,6 +2013,34 @@ mod tests { } } + #[test] + fn rate_test_points_are_the_quartiles_and_midpoint() { + let p = rate_test_points(2.0, 6.0); + assert_eq!(p, [3.0, 4.0, 5.0]); + // A quartile is bitwise the midpoint of the half a bisection makes. + let (a, b) = (0.1, 0.7); + let p = rate_test_points(a, b); + assert_eq!(p[0], rate_test_points(a, p[1])[1]); + assert_eq!(p[2], rate_test_points(p[1], b)[1]); + } + + #[test] + fn rate_error_peak_is_the_vertex_of_the_parabola_at_the_largest_sample() { + // Samples of 1.69 - (x - 1.3)^2 at the quartiles and midpoint of + // 0..4 (zero at x = 0): the largest is at x = 1, and the parabola + // through it and its neighbours x = 0 and 2 has its vertex at 1.3. + let g = |x: f64| 1.69 - (x - 1.3) * (x - 1.3); + let v = rate_error_peak(0.0, 4.0, [g(1.0), g(2.0), g(3.0)]).unwrap(); + assert!((v - 1.3).abs() < 1e-12, "{v}"); + // The same with the sign reversed (an interpolant below the rate). + let v = rate_error_peak(0.0, 4.0, [-g(1.0), -g(2.0), -g(3.0)]).unwrap(); + assert!((v - 1.3).abs() < 1e-12, "{v}"); + // A symmetric profile peaks at the midpoint; no samples, no peak. + let v = rate_error_peak(0.0, 4.0, [0.75, 1.0, 0.75]).unwrap(); + assert!((v - 2.0).abs() < 1e-12, "{v}"); + assert_eq!(rate_error_peak(0.0, 4.0, [0.0; 3]), None); + } + #[test] fn a_window_above_the_grid_adds_no_rows() { let model = SteepRate(steep);