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);