Simulation Engine

The expensive step of a Bloch wave simulation is the eigen-decomposition of the structure matrix A, one matrix per incident beam direction (tilt). PyExtal parses the .dat file, computes structure factors, builds the tilt grid and selects the strongly excited beams in its compiled Bloch engine (pyextal.blochwave), and solves the eigenproblems with NumPy (pyextal.blochNumpy) on a thread pool.

Every simulation entry point takes a backend argument. backend='numpy' is the default and the only backend in the published wheels, so you never need to pass it.

import pyextal.callBloch as cb

param = cb.bloch('Si_systematic.dat', ncores=8)
I = cb.LARBED_tilt(param, 3000.0, param.nout)

CBEDROI, CoarseOptimize and FineOptimize take a backend argument too and pass it to every simulation they run.

Number of cores

ncores sets the number of worker threads.

  • cb.bloch, cb.bloch_run, cb.tilt_run and cb.simulate default to ncores=1. Pass ncores explicitly when calling them yourself.

  • The refinement classes use constant['NCORES'] from pyextal.Constants (default 8).

  • ncores=0 means “all physical cores”.

How the eigensolver works

bloch_run_numpy / tilt_run_numpy in pyextal.blochNumpy run these stages:

  1. Parse (bloch_parse, compiled engine): crystal, structure factors Ugh, beam list.

  2. Tilt grid (_generate_tilt_grid): the incident directions inside the aperture set in the .dat file.

  3. Beam selection (_build_all_A_matrices, pass 1): excitation errors and gn factors are vectorised over all tilts. The strongly excited beams of each tilt are chosen by the reentrant engine routine beamcheck_batch, called on the thread pool.

  4. A-matrix build (_build_all_A_matrices, pass 2): tilts with the same ndiag are processed in chunks. Each tilt’s Ugh sub-blocks are gathered in one batched indexing call, and the second-order Bethe correction is one batched matmul per chunk. Everything is float32/complex64.

  5. Eigensolve (_numpy_eigensolve): numpy.linalg.eig on each exact-size (ndiag, ndiag) matrix, with no padding.

  6. Assembly (_assemble_simparams): eigenvectors are scaled by gn, and C⁻¹ is computed with batched numpy.linalg.inv. Results are stored as Data objects in param.tilts.

Stages 3 to 6 run on a ThreadPoolExecutor. NumPy releases the GIL inside LAPACK, matmul and large array operations, so plain threads run in parallel without the pickling cost of multiprocessing.

BLAS threading

Each worker thread must use one BLAS thread. Otherwise every LAPACK call inside the pool starts its own team of OpenBLAS or MKL threads, and the CPU is oversubscribed. On a 12-core machine this made an 8-thread pool no faster than serial, and a 12-thread pool slower than serial. PyExtal pins BLAS to one thread around the pool with threadpoolctl, a required dependency. If threadpoolctl is missing, the engine falls back to mkl-service. If neither is available it emits a RuntimeWarning and runs unpinned.

Performance

Measured on an Intel i9-12900K (12 physical cores) under WSL2, with OpenBLAS as NumPy’s BLAS. Times are the best of three runs of a full bloch() call (parse, beam selection, eigensolve, assembly).

At ncores=8:

Case

Mean ndiag

Tilts

Time

CBED Si systematic

166

1681

4.55 s

LARBED Si [110]

69

2821

1.36 s

LARBED Cu₂O

89

1961

2.58 s

LARBED Si [133]

40

5025

1.18 s

LARBED GGG

136

441

1.48 s

The refinement loop (tilt_run followed by LARBED_tilt) scales the same way.

Scaling with core count, full bloch():

Cores

CBED Si

LARBED Si [110]

1

30.5 s

7.56 s

2

15.9 s

4.01 s

4

8.20 s

2.18 s

8

4.55 s

1.36 s

12

3.93 s

1.24 s

Where the time goes for CBED Si at 8 cores: beam selection plus A-matrix build 1.8 s, eigensolve 2.5 s, assembly 0.4 s, parsing under 0.1 s.

Inputs with thousands of small matrices (LARBED Si [133]: 5025 tilts, ndiag ≈ 40) are dominated by Python overhead per tilt rather than LAPACK work, so they gain least from more cores.