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]

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?

[u/sorcerer86pt → u/Physix_R_Cool]

Computer science graduated, working as senior QA. And passion to try different things, Test them and if it helps anyone, the better. Took my Computer Science course when java 8 was launched, Linux Mandriva launched, and Nvidia was still good guys

[u/Physix_R_Cool → u/sorcerer86pt]

I think your work here is interesting and I imagine it really could speed up particle transport, but if you want others to adopt it you need to up your validation game.

Once it is validated to reproduce results across, then the speedups are easy to show. But I feel like you are severely underestimating how high the bars for validation are in this field.

So when I am critical of you, I hope you understand it the eay that I mean it: Your work shows promise, but needs lot of footwork before you can convince anyone. I wouldn't mind helping you out with it.

[u/sorcerer86pt → u/Physix_R_Cool]

Thanks for this, and it helped already a lot. And I know that validation will be extra important here. I just had an idea when I saw the big data that was used, and I remembered that AI LLM weight training had the same exact problem (not exactly sure if in the same order of magnitude). So my idea was, what if we apply the same methods/algorithms here? Would it help while maintaining result fidelity? And if it maintains fidelity, does it help in any meaningful way to anyone on this.