# 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.

```python
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`](https://github.com/joblib/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.