4163 inconsistent bosch hale correction in beam fusion calculation by grmtrkngtn · Pull Request #4164 · ukaea/PROCESS

grmtrkngtn · 2026-04-07T15:37:05Z

Description

The beam fusion model is used to calculation fusion power from beam/plasma interactions.

beam_fusion() -> beamcalc() -> alpha_power_beam() -> bosch_hale_reactivity()

bosch_hale_reactivity should be passed a profile, currently it is passed a volume averaged temperature.

This results in a much higher value of beam_fusion being calculated than might be expected and is inconsistent with the way the function should work according to its own doc-string and reference..

This is also inconsistent with how the calculation is done elsewhere in the code, where fusion_rate_integral() is used.

This affects a correction ratio for the reactivity in alpha_power_beam() and has a significant impact on the beam fusion power calculated by PROCESS.

To amend this, we should pass in an ion temperature profile directly; and integrate this as is done elsewhere in the code, and use this integral to calculate the correction ratio.

With this change PROCESS now produces a beam fusion power in line with another code - METIS - for a machine of similar size and operation.

Checklist

I confirm that I have completed the following checks:

My changes follow the PROCESS style guide
I have justified any large differences in the regression tests caused by this pull request in the comments.
I have added new tests where appropriate for the changes I have made.
If I have had to change any existing unit or integration tests, I have justified this change in the pull request comments.
If I have made documentation changes, I have checked they render correctly.
I have added documentation for my change, if appropriate.

Pass plasma profiles to the beam fusion power calculation to ensure consistent reactivity evaluation.

…cstrings. Passed plasma_profile. Renamed and refactored beamcalc()

…ocstring terms

codecov-commenter · 2026-04-08T16:02:10Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 49.55%. Comparing base (09c52f5) to head (4b7d80e).
⚠️ Report is 10 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #4164      +/-   ##
==========================================
+ Coverage   48.04%   49.55%   +1.51%     
==========================================
  Files         144      148       +4     
  Lines       30211    29736     -475     
==========================================
+ Hits        14515    14736     +221     
+ Misses      15696    15000     -696

☔ View full report in Codecov by Sentry.
📢 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.

chris-ashe · 2026-04-09T13:59:24Z

process/models/physics/fusion_reactions.py


-        - J. W. Sheffield, “The physics of magnetic fusion reactors,” vol. 66, no. 3, pp. 1015-1103,
-        Jul. 1994, doi: https://doi.org/10.1103/revmodphys.66.1015.
+        - J. W. Sheffield, “The physics of magnetic fusion reactors,”


Can I ask that the clickable DOI link is preserved

chris-ashe · 2026-04-09T14:13:46Z

process/models/physics/fusion_reactions.py

+        fusion_rate_integral(
+            plasma_profile,
+            BoschHaleConstants(**REACTION_CONSTANTS_DT),
+        ),


I am confused as to how this is different from what was already being done. Also have you weighted the densities and temperatures to match those for the ions and not the electrons?

# Calculate the fusion reaction rate integral using Simpson's rule sigmav = integrate.simpson( fusion_rate_integral(self.plasma_profile, dt), x=self.plasma_profile.neprofile.profile_x, dx=self.plasma_profile.neprofile.profile_dx, )

We aren't changing the assumption about the ion temperature and density profiles. In the old version and the new version they are both scaled by the volume averaged quantities (is there a better way to do this?)

The real change of this PR is:

ratio = ( sigmav_dt / bosch_hale_reactivity( np.array([temp_plasma_ion_vol_avg_kev]), BoschHaleConstants(**REACTION_CONSTANTS_DT), ).item() )

We were passing volume average temp instead of the profile.

chris-ashe · 2026-04-09T14:14:12Z

process/models/physics/fusion_reactions.py

+    # and the Bosch-Hale reactivity evaluated with the same profile integral.
+    sigmav_dt_correction = sigmav_dt_average / bh_profile_average
+
+    hot_deuterium_rate = 1e-4 * beam_reaction_rate(


Why 1e-4 in these cases?

I am not sure to be honest, this is from the old code not something I introduced. I'll take a look into it.

I believe it's to convert from cm^{2} to m^{2}.

_beam_fusion_cross_section outputs in cm^{2}

Ah ok, Can I ask that the variable name is changed to represent what it really is. Looking at the docstring of beam_reaction_rate it says it gives out float: Hot beam fusion reaction rate (m^3/s).. In this case the fusden_ prefix should be used. Also if it actually in per cm cubed then please add _cm to the end of the name

I've moved this conversion into the beam_reaction_rate function. I think this is better as it now outputs SI units. Also updated the docstring.

chris-ashe · 2026-04-09T14:15:35Z

process/models/physics/fusion_reactions.py

-        range of beam velocities up to the critical velocity.
-
-    References:
-        - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and thermal reactivities,”


Why have references been removed

Mistake, I've added them back in.

timothy-nunn · 2026-04-09T16:01:54Z

Once @chris-ashe approves I will check this out and do a code review

chris-ashe

Looks good logically. Are we not adding the profile correction for the beam width? Should be easy enough to do.

The docs also need to be updated.

I am not asking to replace everything, but can new variables you add please follow the style guide

chris-ashe · 2026-04-10T09:44:45Z

process/models/physics/fusion_reactions.py

+    # and the Bosch-Hale reactivity evaluated with the same profile integral.
+    sigmav_dt_correction = sigmav_dt_average / bh_profile_average
+
+    hot_deuterium_rate = 1e-4 * beam_reaction_rate(


Ah ok, Can I ask that the variable name is changed to represent what it really is. Looking at the docstring of beam_reaction_rate it says it gives out float: Hot beam fusion reaction rate (m^3/s).. In this case the fusden_ prefix should be used. Also if it actually in per cm cubed then please add _cm to the end of the name

chris-ashe · 2026-04-10T09:45:53Z

process/models/physics/fusion_reactions.py

+        constants.M_TRITON_AMU, tritium_critical_energy_speed, e_beam_kev
+    )
+
+    deuterium_beam_alpha_power = alpha_power_beam(


As per style guide. p_beam_deuterium_dt is the recommended name`

i changed hot_deuterium_rate to sigv_beam_deuterium. It's not a fusion density so I don't think fusden is right.

I changed quite a few additional variable names

process/models/physics/fusion_reactions.py

grmtrkngtn · 2026-04-10T15:06:46Z

Looks good logically. Are we not adding the profile correction for the beam width? Should be easy enough to do.

I would prefer to keep this PR for the correction + refactor, then make a separate issue for the beam width (new phyiscs)

… commented. updated function call and tests. renamed hot_deuterium_rate to sigv_beam_deuterium

…quent functions

grmtrkngtn · 2026-04-10T16:01:10Z

The docs also need to be updated.

Updated the docs now

grmtrkngtn added 3 commits April 7, 2026 14:04

Fix Bosch–Hale reactivity calculation in alpha_power_beam

58108d4

Pass plasma profiles to the beam fusion power calculation to ensure consistent reactivity evaluation.

update alpha_power_beam docstring

6dc56f1

incorporate changes to beam_fusion calculation in physics.py

cbb4576

grmtrkngtn requested a review from a team as a code owner April 7, 2026 15:37

grmtrkngtn linked an issue Apr 7, 2026 that may be closed by this pull request

Inconsistent Bosch-Hale correction in beam fusion calculation #4163

Open

chris-ashe self-assigned this Apr 8, 2026

grmtrkngtn added 4 commits April 8, 2026 15:58

Re-work beam_fusion execute routine in a more modular way. Updated do…

62ff5da

…cstrings. Passed plasma_profile. Renamed and refactored beamcalc()

updated alpha_power_beam() function with new refactor. updated some d…

9cbf230

…ocstring terms

update fusion reaction tests to use refactored functions

91b4c16

fixed inputs to beam_fusion, plasma profile in the wrong place

e85d2a4

grmtrkngtn added 2 commits April 9, 2026 10:35

removed debugging print

d67dc6d

removed debugging line

f683a23

timothy-nunn unassigned chris-ashe Apr 9, 2026

timothy-nunn requested review from chris-ashe April 9, 2026 12:57

timothy-nunn self-assigned this Apr 9, 2026

chris-ashe reviewed Apr 9, 2026

View reviewed changes

grmtrkngtn added 2 commits April 9, 2026 15:34

Added references to beam_slowing_down_state

b22a368

Added back clickable references to beam_fusion

616b973

chris-ashe requested changes Apr 10, 2026

View reviewed changes

grmtrkngtn added 2 commits April 10, 2026 15:56

updated deuterium beam power variable name

c530e45

updated tritium beam power variable name

db7ee5c

grmtrkngtn added 5 commits April 10, 2026 16:09

beam_reaction_rate docstring to correct otput units

7154405

Updated beam_reaction_rate, moved unit conversion inside function and…

711df66

… commented. updated function call and tests. renamed hot_deuterium_rate to sigv_beam_deuterium

renamed sigmav_dt_correction to f_sigmav_dt_profile_correction

d6b0eb5

renamed bh_profile_average to sigv_dt_bosch_hale_profile_avg

089a691

Renamed beam_slow_time to t_beam_slow

74f8da1

grmtrkngtn added 8 commits April 10, 2026 16:27

renamed deuterium_density

45c1ede

renamed tritium_density

44c17f1

renamed hot_beam_density to nd_beam_hot

59d94fb

renamed beam_depo_energy to e_beam_deposited_kev

cee64c9

renamed alpha_power_beam outputs

1993de3

update call of beam_fusion

3bde93c

fixed beam fusion call in stellarator

34c3c20

update beam_reactions doc to reflect changes in beam_fusion and subse…

4b7d80e

…quent functions

Conversation

grmtrkngtn commented Apr 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Checklist

Uh oh!

codecov-commenter commented Apr 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

timothy-nunn commented Apr 9, 2026

Uh oh!

chris-ashe left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

grmtrkngtn commented Apr 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

grmtrkngtn commented Apr 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

grmtrkngtn commented Apr 7, 2026 •

edited

Loading

codecov-commenter commented Apr 8, 2026 •

edited

Loading

grmtrkngtn commented Apr 10, 2026 •

edited

Loading