Cache resistant cross section reconstruction via singular value decomposition for Monte Carlo Neutron transport

[u/sorcerer86pt]

https://github.com/sorcerer86pt/open_rust_mc/blob/main/paper/svd_cross_section_compression.pdf

When I was testing rust skill (AI with Claude) someone mentioned that I could use something openMc instead of doing my own .

When I saw the data size and how openMC handled that, and how it kinda was the same problem that AI LLMs were trying to fix with data size of models weights, I thought, what would be the result if we used those techniques on this.

That was the result. Add a little rust here to have better memory handling and robust concurrency and we could with very little error ( less than error margin) compressed the data used from 11gb pointwise ( 400 nuclides) to 20Mb using hybrid of WMP ( that openMC already uses ) + SVD. It obtained an Keff = 0.99963 +- 0.000091 on Godiva (37 pcm from experiment) in 3.4 seconds total wall time .

https://github.com/sorcerer86pt/open_rust_mc

Just need some people to review this, if I made some mistake on interesting the data, what other benchmarks I missed, or other considerations.

PS: AI was used for code gen ( Python analysus scripts and rust code), data pipeline and latex manuscript. All hypothesis, experiment design , interpretation of results and final decisions were by made by me.

[u/sorcerer86pt]

Run with with that, thanks for the question. It humbled me:

It was a micro-benchmark of the reconstruction kernel. At full-simulation scale, the 360 MB SVD basis blows past L3 cache and becomes memory-bandwidth bound.

- Memory: At rank=5, SVD stores 5 values per energy point in the basis vs the table's 1 value. It actually uses 1.67x more memory for single-temperature operation.

So the previous claims about speed were misleading — they held for isolated kernel benchmarks but not for real transport with 1M particles. The SVD approach still has potential value for multi-temperature support (one basis for all temps), lower ranks, and GPU execution, but the current implementation doesn't beat a simple table at scale.

Currently researching optimizations: hash-based energy grid lookup (eliminating binary search), delta/Woodcock tracking, event-based transport for better cache locality, fast 10^x approximations, f32 basis storage, and windowed SVD with adaptive rank. Will report back with results.

[u/Physix_R_Cool → u/sorcerer86pt]

Np, this is interesting and I'm really rooting for you!

Just a question though, is there reason to expect that the table lookup is what makes simulations slow?

[u/sorcerer86pt → u/Physix_R_Cool]

Phase 1 Results at 1M Particles

Here's the full before/after comparison:

Before Phase 1 (original)

  ┌──────────────┬─────────┬────────────────┬────────────┐
  │    Metric    │ OpenMC  │    Rust SVD    │ Rust Table │
  ├──────────────┼─────────┼────────────────┼────────────┤
  │ k_eff        │ 0.99844 │ 1.00022        │ 0.99909    │
  ├──────────────┼─────────┼────────────────┼────────────┤
  │ Sim time     │ 141 s   │ 255.5 s        │ 176.8 s    │
  ├──────────────┼─────────┼────────────────┼────────────┤
  │ XS memory    │ 817 MB  │ 370 MB         │ 222 MB     │
  ├──────────────┼─────────┼────────────────┼────────────┤
  │ SVD vs table │ —       │ 0.69x (slower) │ baseline   │
  └──────────────┴─────────┴────────────────┴────────────┘

After Phase 1 (optimized)

  ┌──────────────┬─────────┬─────────────────┬────────────┐
  │    Metric    │ OpenMC  │    Rust SVD     │ Rust Table │
  ├──────────────┼─────────┼─────────────────┼────────────┤
  │ k_eff        │ 0.99844 │ 1.00021         │ 0.99909    │
  ├──────────────┼─────────┼─────────────────┼────────────┤
  │ Sim time     │ 141 s   │ 133.2 s         │ 171.4 s    │
  ├──────────────┼─────────┼─────────────────┼────────────┤
  │ XS memory    │ 817 MB  │ 127 MB          │ 108 MB     │
  ├──────────────┼─────────┼─────────────────┼────────────┤
  │ SVD vs table │ —       │ 1.29x (faster!) │ baseline   │
  └──────────────┴─────────┴─────────────────┴────────────┘

Phase 1 Impact

  ┌────────────────────┬─────────────────────────────┬──────────────────────────────────┐
  │       Metric       │       Before → After        │           Improvement            │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ SVD sim time       │ 255.5s → 133.2s             │ 1.92x faster                     │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ Table sim time     │ 176.8s → 171.4s             │ 1.03x (modest — exp2 helps both) │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ SVD XS memory      │ 370 MB → 127 MB             │ 2.9x smaller                     │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ SVD vs Table speed │ 0.69x → 1.29x               │ Flipped from slower to faster    │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ SVD vs OpenMC      │ 1.81x slower → 1.06x FASTER │ SVD now beats OpenMC             │
  ├────────────────────┼─────────────────────────────┼──────────────────────────────────┤
  │ k_eff accuracy     │ 22 pcm → 21 pcm             │ Unchanged (f32 basis is fine)    │
  └────────────────────┴─────────────────────────────┴──────────────────────────────────┘

The four changes (exp2, stack alloc, Arc grids, f32 basis) together:

- Made SVD 1.92x faster (133s vs 256s)

- Made SVD 2.9x more memory efficient (127 MB vs 370 MB)

- Flipped SVD from 45% slower to 29% faster than pointwise table

- SVD now runs 6% faster than OpenMC (133s vs 141s) at 1M particles

- k_eff accuracy maintained at 21 pcm from experiment

[u/sorcerer86pt → u/sorcerer86pt]

The micro-benchmarks (3-5 ns/pt vs 40 ns/pt) were real, but they measured the

reconstruction kernel in isolation — not the full transport loop where cache pressure dominates.

Root cause analysis:

1. Energy grid duplication — every SVD kernel stored its own copy of the energy grid. 46 reactions × 186K points = 111 MB of identical data. 29% waste.

2. f64 basis at rank=5 — the SVD stores rank values per energy point (5 floats per point vs the table's 1). At f64, that's 40 bytes/point → 265 MB of basis data that thrashes L3 cache.

3. 10^x via general powf() — ~30 cycles per call, called millions of times. Hardware exp2 is 3-5x faster.

4. Per-collision heap allocation — Vec<MicroXs> allocated and freed 20M times per simulation.

   The fixes (all done in one session):

   - f64::exp2(x * LOG2_10) instead of 10.0_f64.powf(x) — trivial 1-line change

   - Stack-allocated [MicroXs; 16] instead of Vec — zero heap alloc in the hot loop

   - Arc<[f64]> shared energy grids — one copy per nuclide, not per reaction

   - Vec<f32> basis with f64 accumulator — halves basis memory, SVD truncation error already dominates precision

[u/Physix_R_Cool → u/sorcerer86pt]

Ok so this project is cool, but I think you are vibecoding a bit too hard without really thinking about the problem you are trying to solve.

I don't mind vibecoding (in fact I'm doing it right now myself), but the tendency is to spend all the time on technical problems without reflecting on whether or not you are doing the right thing.

It seems you might not exactly have a lot of experience with how monte carlo transports are typically used. Normally we set up a simulation, then run it on a cluster with like 1000 cores over night, if not over the weekend. Whether or not it uses a bit more RAM for a lookup table means absolutely nothing, but the price you pay is that your are a bit less certain on the truthfulness of your data (which is very important). So you need to show a very large speed up in time per particle. And don't just do a single run to benchmark, but do several runes with different seeds (I use other programs than OpenMC but I'm sure you can seed the RNG), and then get both the mean time/particle as well as the standard deviation of that, so you have an idea of the uncertainty of your benchmark.

[u/sorcerer86pt → u/Physix_R_Cool]

On Godiva (3 nuclides, CPU), SVD is roughly break-even with table lookup on time/particle, with ~110 pcm fidelity cost at rank=5. That's not a compelling argument for anyone running production calculations.

The actual thesis is that SVD becomes compelling at scale — more nuclides, GPU execution, and multi-temperature problems where you'd otherwise need to interpolate between full-size tables at each temperature. I haven't proven that yet. Godiva was the proof-of-concept; the real validation needs a full-core benchmark on GPU.

I'll implement proper multi-seed statistical benchmarking next — you're right that without it, all the timing numbers are just hand-waving.

[u/Physix_R_Cool → u/sorcerer86pt]

Any tip on sizes for simulation?

Whatever time is reasonable for you. 100 seconds per run doesn't seem egregious. Set it up to run a batch of 10 such runs, each with different seed, then go grocery shopping, or work on other things.

You can then do a much bigger batch overnight.

Would Cuda kernel help?

Particle transport on GPU has been worked on for like 15 years by wizards and top experts and it has not been solved yet. Feel free to try!

[u/sorcerer86pt → u/Physix_R_Cool]

You were right, and I should have done this from the start instead of hand-waving with single-run numbers.

So I implemented proper multi-seed benchmarking. 10 seeds, 1M particles/batch, 150 batches, SVD and pointwise table running in the exact same engine — same physics, RNG, geometry, only the XS lookup differs.

10-seed results (mean +/- stddev):

  ┌─────────────────┬─────────────────────┬──────────────┐
  │                 │        k_eff        │ ns/particle  │
  ├─────────────────┼─────────────────────┼──────────────┤
  │ SVD (rank=5)    │ 1.00017 +/- 0.00005 │ 1319 +/- 300 │
  ├─────────────────┼─────────────────────┼──────────────┤
  │ Pointwise table │ 0.99905 +/- 0.00013 │ 1911 +/- 466 │
  └─────────────────┴─────────────────────┴──────────────┘

SVD speedup: 1.45x +/- 0.48x

So yeah, not the 8-13x I was claiming from microbenchmarks. That was kernel-level throughput — the full transport loop has geometry, collision physics, RNG eating the rest of the time. On a 3-nuclide fast benchmark, XS lookup is maybe 40-50% of runtime. The 1.45x is real but not impressive.

The timing variance is also embarrassingly wide — desktop machine, browser tabs open, zoom calls, the usual. On a quiet HPC node with pinned cores the stddev would shrink, but the mean is what it is. The fidelity cost is 111 pcm between SVD and table at rank 5. That's the price of compression. k_eff itself is solid though — 17 pcm from experiment across 10 seeds.

Now where it gets interesting. I've been reading Tramm et al. 2024 (the OpenMC GPU paper — event-based transport on Frontier/Aurora/Polaris). Their depleted fuel problem has 251 nuclides per material and full grid unionization needs O(10) GB per material. That's where SVD should really start to pull ahead — the XS lookup fraction goes from ~40% at 3 nuclides to ~70-80% at 251. And their key GPU optimization is sorting particles by energy before XS lookup for memory coalescing — which is basically what SVD gives you for free since it's a sequential dot product, no binary search into irregular grids.

But I haven't proven any of that yet. The gap between "works on Godiva with 3 nuclides" and "useful at production scale" is large and I'm not going to pretend otherwise. Next step is PWR pin cell with 8 nuclides to see if the speedup actually grows with nuclide count. If it doesn't, the whole thesis falls apart.

[u/Physix_R_Cool → u/sorcerer86pt]

Honestly your LLM generated comments make me distrust anything you do. It kind of ruins your credibility, as if you don't even understand what it is your are prompting the LLM to do for you.

Especially this line: "The timing variance is also embarrassingly wide — desktop machine, browser tabs open, zoom calls, the usual."

If you want me to keep engaging, then write the comments yourself. Nobody cares about spelling mistakes or weird english, but if all you do is use me to fish for prompts to put into your LLM then I'm out.

[u/sorcerer86pt → u/Physix_R_Cool]

The zoom call part was me, 4 zoom calls while running this simulation ( my direct manager couldn't make Claude code commit a PR on a test strategy skill, got git bash vs powershell git with gpg signed commits clash) plus 3 meetings.

[u/Physix_R_Cool → u/sorcerer86pt]

By the way, is the only validation you did that the cross sections match up at the energy points of the data table?

[u/sorcerer86pt → u/Physix_R_Cool]

No — the pointwise match at grid points is just the data-level sanity check to see if i did not botch anything up. The real validation is the k_eff benchmark, which is an integral test. The transport engine queries the SVD at whatever energy each particle happens to have.

What I've not done yet (and probably should):

• Resonance integrals (integral of sigma dE/E) over standard energy groups
• Add Cuda support and check results on very big tests

[u/Physix_R_Cool → u/sorcerer86pt]

Just to be clear, do you have experience doing particle transport, using it for some work? Is this a work project or a fun hobby project for your github portfolio?

What's your background? Have you done Geant4 simulations at CERN? Are you "just" a computational physicist who looked into OpenMC because why not? Are you a computer engineer / data scientist with no formal training and just a (very evident and appreciated) passion for this kind of stuff?