Optimise neutral_mixed flux limiter implementation

[mikekryjak]

Purpose

I have been trying to merge the separated flux limiters for several years (see #556, #451, mikekryjak#15 and mikekryjak#6). My efforts were hampered by extreme performance degradation which I couldn't understand.

I have recently been conferring with LLMs to try and understand how a code could slow down by 100x just because we limit conduction separately. I found that the code is extremely sensitive to regularisation issues (e.g. denominators going to zero, non-differentiable things, etc). There a lot of small things in our existing implementation and my implementation of the separate limiters which caused problems.

I decided to separate the PRs and implement the regularisation first as well as the new limiter features apart from parallel and split limiters (so smooth SOLPS like limiter form and new Vth definition).

The new changes give robustness improvements a scan of flux_limiter values, while being much closer to the unregularised results. They will change the physics result to be more accurate. Previous profiles can be recovered by tuning the flux limiter.

Change Summary

• Thermal speed definition has corrected, which results in an offset to flux_limit. The new default of 0.5 is equivalent to the old default of 0.2. If you have simulations where flux_limit is not set to default, your results will be affected.
• neutral_lmax is now user-settable and only controls the collisionality floor for separation of concern.
• diffusion_limit now caps Dmax even if flux_limit is used, as opposed to replacing the flux limitation with a D limitation.
• Collisionality floor now an actual softFloor instead of addition. Nu is no longer always higher by Rnn everywhere and there is no longer a kink at the floor.
• Rnn is now called nu_pseudo_mfp.
• Changed default collisions mode to AFN as the more accurate model.
• Implemented SOLPS-style variable sharpness limiter (currently defaults to 1, consistent with previous behaviour)
• Dmax denominator now correctly uses the perpendicular instead of total gradient.
• abs(Grad_perp(logPnlim)) in Dmax denominator creates a non-differentiable kink. Replaced with smooth sqrt(abs(Grad_perp(logPnlim))^2 + gradient_floor^2) and removed the now redundant 1/neutral_lmax` floor.
• Add a smooth Grad_perp(logPnlim) ceiling to prevent crashes in steep transients when the limiter is saturated.
• Many many runs on my fork (https://github.com/mikekryjak/hermes-3/tree/neutral-fluxlim-full)
• Final validation runs of clean PR
• Documentation
• Reran and updated 2D-production and 2D-recycling tests
• Added new unit tests for the flux limiter

Validation

I've ran comprehensive test on my low-res ST40 case and performed a second verification on a MAST-U case. I did a scan in each of the normalisation parameters to figure out sensible defaults.

Results of the runs:

• neutral_lmax set to physical 1m. It now only floors collisionality. Increasing it doesn't change the answer much.
• limiter_gradient_ceiling is at 0.1: only mildly activated in steady state, but prevents crashes at initial transients when the flux is heavily limited.
• limiter_gradient_floor is at 0.01: prevented huge slowdowns due to shallow gradient regions when separate limiters are used (next PR).
• Changing Grad() to Grad_perp had no significant impact but is more correct.

My low res case scan

Profiles:

Performance:

I found that a lot of the speed variability I've been getting is due to the slope limiter choice. Here are the scans repeated with the Van Albada slope limiter which shows that the new changes now have very consistent performance and are still better than master:

AI Assistance

• I wrote all the code, comments and text in this PR myself.
• Used LLMs extensively in investigating and planning the changes.
• Used LLMs extensively to check my code and simulation inputs for typos, bugs and poor C++ practice.
• Sometimes used LLMs to write code on a separate branch for prototyping.

Documentation

The neutral_mixed section was rewritten to include the new changes.

Review Notes

• This PR will change the result. Firstly, the new Vth definition means that e.g. a flux limiter value of 0.5 is equivalent to the old 0.2. I have kept the defaults to be equivalent to old behaviour. Secondly, the regularisation is now less impactful on the simulation, so the profiles will change, generally towards less limitation. This can be compensated for with a different flux limiter setting.
• Appreciate comments on setting/variable names.

[codecov]

Codecov Report

❌ Patch coverage is 82.35294% with 9 lines in your changes missing coverage. Please review.
✅ Project coverage is 58.61%. Comparing base (731e644) to head (5dcbdc8).
⚠️ Report is 25 commits behind head on master.

Files with missing lines	Patch %	Lines
src/neutral_mixed.cxx	82.35%	7 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #605      +/-   ##
==========================================
+ Coverage   56.92%   58.61%   +1.68%     
==========================================
  Files          97       98       +1     
  Lines       10092    10274     +182     
  Branches     1462     1471       +9     
==========================================
+ Hits         5745     6022     +277     
+ Misses       3757     3640     -117     
- Partials      590      612      +22     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mikekryjak]

I've been using the prek pre-commit hook and fixing formatting as I go. Things passed locally, not sure how clang-format is still complaining...

[malamast]

Hi @mikekryjak , I'll try to test this in the next few days but do you think we should split the flux limiters into parallel and perp? I don't have any strong opinion.

[mikekryjak]

Hi @mikekryjak , I'll try to test this in the next few days but do you think we should split the flux limiters into parallel and perp? I don't have any strong opinion.

Yes, we should, but I wanted that in a separate PR. This one keeps the same limiter form as master.

[mikekryjak]

I realised that leaving limiter_gradient_ceiling and limiter_gradient_floor in normalised units is pretty bad, because they will then depend on the normalisation parameters. I just made them normalised (units of 1/m).