Compare grid resolution and dynamical backends

Compute nickel master patterns at two grid resolutions, project onto an identical detector, and quantify the interpolation difference. Then compare the scattering-matrix backend with the Bloch eigensolver at identical settings.

An explicitly prescribed exponential depth distribution keeps this numerical comparison short and deterministic. It is an approximate depth model, not a replacement for the native Monte Carlo workflow in the preceding examples. Two grids demonstrate a convergence procedure; agreement between two solvers is not evidence that the physical model is converged.

from time import perf_counter

import anyplotlib as apl
import numpy as np
import torch
from orix.quaternion import Rotation

import de_ebsd as de

Hold the experiment and depth model fixed

Use double precision in both solvers for this comparison. Changing only npx preserves the crystal, reflection cutoff and structure factors.

device = "cpu"
crystal = de.Crystal.from_parameters(
    a=0.35236,
    space_group=225,
    atoms=[("Ni", (0.0, 0.0, 0.0))],
    debye_waller=0.0035,
)
microscope = de.Microscope(energy=20e3, sample_tilt=70.0, probe_diameter_nm=10.0)
detector = de.EBSDDetector(
    shape=(80, 100),
    pattern_center=(0.45, 0.65, 0.60),
    sample_tilt=microscope.sample_tilt,
    tilt=5.0,
    px_size=400.0,
)
depth = de.DepthDistribution.single_energy(
    microscope.energy_kev,
    thickness=40.0,
    depth_step=2.0,
    decay_length=15.0,
)

Inspect the crystal

This is the actual expanded Ni unit cell in Cartesian nanometres. Drag to rotate the native anyplotlib figure; the same figure is usable in SpyDE.

crystal.show()

Inspect the microscope, probe and detector

The beam is vertical in the laboratory side view. Positive sample tilt slopes the surface down toward the detector; the arrow n points out of the sample. Incidence is measured from n (70 degrees), not from the surface (20 degrees). The example’s assumed 10 nm circular parallel probe projects to a 29.24 x 10 nm footprint; it does not model electron scattering or broadening.

Here the assumed effective screen pitch is 400 micrometres per binned pixel: the 100 x 80 image spans a 40 x 32 mm screen, with a 19.2 mm normal distance to its pattern center. This is screen calibration, not camera-sensor pitch. PCy is measured down from the top; PCz is normalized by screen height. The screen is tilted 5 degrees from vertical, 25 degrees relative to the sample plane. Its center is 3.11 mm above the impact point in laboratory height, and 19.43 mm above the sample plane along the sample normal. Click the detector face or edit the numeric fields to explore the geometry.

Edits change this inspector only. In a notebook, explicitly read setup.to_microscope() and setup.to_detector(), rebuild the simulator, and rerun downstream cells to use them. Exported gallery controls cannot recompute Python results. detector.show(microscope=microscope) offers the same geometry inspector when starting from a detector object.

setup = microscope.show(detector=detector, crystal=crystal)
setup

Inspect the depth slices used by both backends

The prescribed model has 20 depth bins, each 2 nm thick, spanning 0–40 nm below the surface. The viewer reads the real distribution and the same 99% cutoff used by backend preparation. Select a slice to inspect its yield and cumulative fraction; selection never changes the distribution.

These are normal-depth integration slices, not crystal layers or a simulated electron trajectory. The stored exponential is sampled at z = j * dz; PRISM evaluates the propagated wave at (j + 1) * dz. Both solvers use the absorption-corrected weights internally; this viewer shows uncorrected exit statistics, so it does not mistake an increasing correction for deeper yield.

depth.show(energy_kev=microscope.energy_kev)

Prepare the same diffraction model at two grid sizes

coarse_sim = de.MasterPatternSimulator(
    crystal,
    microscope,
    npx=100,
    dmin=0.05,
    monte_carlo=depth,
    backend="bloch_prism",
)
coarse_sim.structure_factors()
fine_sim = coarse_sim.refine(npx=200)
assert fine_sim.structure_factors() is coarse_sim.structure_factors()

Compute both grids

These timings include preparation and any first-use overhead, so they are illustrative rather than a controlled solver benchmark.

start = perf_counter()
coarse = coarse_sim.compute(device=device, dtype=torch.float64, chunk_size=512)
print(f"201 x 201 per hemisphere: {perf_counter() - start:.3f} s")
start = perf_counter()
fine = fine_sim.compute(device=device, dtype=torch.float64, chunk_size=512)
print(f"401 x 401 per hemisphere: {perf_counter() - start:.3f} s")
fine.show(projection="sphere", texture_height=128, n_grid=64)

201 x 201 per hemisphere: 1.816 s
401 x 401 per hemisphere: 7.671 s

Compare on the same detector grid

Compare raw intensities so independent image normalization cannot hide an amplitude change. Report relative RMS error against the finer projection. A larger Lambert grid samples the same diffraction model more finely; lowering dmin is a separate convergence axis: it adds reflections rather than sampling the existing ones more finely. Both grids here share one dmin=0.05 reflection list, so this comparison isolates sampling alone.

orientation = Rotation.from_euler([[35.0, 45.0, 10.0]], degrees=True)
coarse_image = coarse.project(detector, orientation, rescale=False)[0]
fine_image = fine.project(detector, orientation, rescale=False)[0]
difference = fine_image - coarse_image
relative_rms = difference.square().mean().sqrt() / fine_image.square().mean().sqrt()
assert torch.isfinite(relative_rms)
print(f"Grid-change relative RMS difference: {relative_rms.item():.3%}")
figure, axes = apl.subplots(1, 3, figsize=(900, 300))
for ax, values, title in zip(
    np.ravel(axes),
    (coarse_image, fine_image, difference),
    ("npx = 100", "npx = 200", "Fine minus coarse (own contrast)"),
    strict=True,
):
    panel = ax.imshow(values.detach().cpu().numpy(), cmap="gray", origin="upper")
    panel.set_title(title)
figure

Grid-change relative RMS difference: 15.876%

Compare the two solvers at identical sampling and depth settings

Use the coarser of the two grids: the exact eigensolver costs ~4x the scattering-matrix backend here, and ~20x at npx=200. Reuse the same supplied depth distribution and float64 directions. This checks implementation consistency; neither solver is being compared with measured data here.

reference_sim = coarse_sim.refine(backend="bloch")
reference = reference_sim.compute(device=device, dtype=torch.float64, chunk_size=512)
scale = reference.data.abs().max().clamp_min(1e-30)
relative_max = (coarse.data - reference.data).abs().max() / scale
assert torch.isfinite(relative_max) and relative_max < 1e-7
print(f"PRISM versus Bloch maximum relative difference: {relative_max.item():.3e}")

PRISM versus Bloch maximum relative difference: 2.218e-13

Continue a convergence study

Increase npx again and repeat the detector-space comparison. Separately lower dmin, increase the depth window, and increase the Monte Carlo electron count until the relevant measured quantity is stable. Changing texture_height in .show() only changes the rendering and cannot improve the computed master pattern.

print("Next grid size:", fine_sim.refine(npx=400).npx)

Next grid size: 400