Particle balance constraint updated by chris-ashe · Pull Request #4395 · ukaea/PROCESS

chris-ashe · 2026-06-29T12:32:35Z

This pull request introduces a comprehensive plasma fuelling model to the codebase, including new physical variables, input parameters, constraints, and documentation. The main focus is to enable detailed particle balance and burnup fraction calculations for deuterium, tritium, and helium-3 in the plasma, and to ensure consistency through new constraint equations. Several new variables and iteration parameters are added to support these models, and the documentation is expanded to explain the physical basis and equations.

Key changes:

1. Plasma fuelling and particle balance model:

Added a detailed plasma fuelling model and particle balance equations for deuterium, tritium, helium-3, and alpha particles to the documentation (plasma_fuelling.md). This includes explanations of fuelling efficiency, recycling, and burnup fractions.
Added new physics variables in PhysicsData for burnup fractions and detailed fusion reaction rates, enabling more granular tracking of fuel species and reactions.

2. New input and iteration variables:

Introduced new input parameters in input.py for fuelling efficiency, injected fuel rates, and fuelling fractions for each species.
Added corresponding iteration variables for these new inputs, allowing them to be varied during optimization and scans.

3. Particle balance and consistency constraints:

Implemented new constraint equations (93–97) in constraints.py to enforce particle balance for tritium, deuterium, helium-3, alpha particles, and to ensure fuelling fractions sum to unity. These constraints help maintain physical consistency in the plasma model.
Updated numerics to include these new constraints, increasing the total number of constraints and variables.

4. Documentation and navigation updates:

Added the new plasma fuelling documentation page to the navigation in mkdocs.yml, making the new model easily accessible in the documentation site.

5. Miscellaneous improvements:

Added Avogadro's number to the constants module for use in physical calculations.
Improved output reporting for constraint errors in scan.py.

These changes collectively provide a robust framework for modeling plasma fuelling, tracking individual fuel species, and ensuring physical consistency in fusion plasma simulations.

🎨 Output additions

Expanded and added more values for rates to the fusion reaction summary page:

Added a fuelling summary page that shows the contour graphs of the fuelling solutions, along with burnup data:

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.

j-a-foster

Just a couple of notes, otherwise happy with the changes.

j-a-foster · 2026-06-29T14:42:21Z

+        f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.4f} MW\n"
+        f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.4f} MW\n\n"
+        f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m$^3$/sec\n"
+        f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m$^3$/sec\n\n"


Why are these /sec and not /s?

j-a-foster · 2026-06-29T14:55:44Z

+        linewidths=2,
+    )
+
+    # Plot star for mfile values


Is it worth adding a note or legend that explains what the star means in the PDF?

jonmaddock

Please sort out the PR description: I think it's too verbose and misses/buries the key point which is that you're introducing constraints to enforce density equilibrium of individual ion species. Why are you making this PR?
Have you removed the original molflow_plasma_fuelling_required as an output? I might have missed that
I counted 6 new optimisation parameters and 5 new constraints: as per our conversation, can you describe how this might work in solution mode, i.e. when we require a determined system?
Some plots to demonstrate these changes would be useful: for example how the constraints are accommodated with increasing te or ne, for example. How does this change the current large tokamak solution?
Please don't rebase until the PR is approved, so we can track comments
How does this compare to the existing burnup calculation? Has it been removed?
The created docs were excellent
I'd like this PR to include what equations and parameters should be included in optimisation and solution scenarios. If the solution system is under-determined, how useful is it?
I'm not sure about the fuelling composition constraint and total fuelling rate. Would individual species rates reduce the dimensionality?
Should the recycling fraction and fuelling efficiency be optimisation parameters? (I realise we've discussed this, but I think it should be made clear why these can be used to solve the constraints).

Thanks for the fixed-up commits, this was much easier to review.

jonmaddock · 2026-06-30T14:05:06Z

+
+        Notes
+        -----
+        The fusion rate is multiplied by two to convert from nucleus pairs to particles,


Not multiplied by 2 here.

jonmaddock · 2026-06-30T14:39:24Z

+        * data.physics.vol_plasma
+        * data.physics.f_plasma_fuel_helium3
+    ) / (
+        data.physics.t_energy_confinement


Does He3 follow the energy confinement time rather than the tau_alpha / tau_E = 5 relation?

This is a point we should discuss, I dont see why we wouldn't treat it the same as 4He

jonmaddock · 2026-06-30T15:20:52Z

+f_{\text{fuelling,D}} + f_{\text{fuelling,T}} + f_{\text{fuelling,3He}} = 1.0
+$$
+
+**It is recommended to have this constraint on as it is a plasma consistency model**


I think you need to be more explicit about the system of equations (i.e. all of the above constraints) and the solution parameters (i.e. optimisation parameters) used to solve them. What's should the user do to enforce all of these constraints in their optimisation problem?

chris-ashe · 2026-07-01T07:44:11Z

@timothy-nunn If I try and use the functions in fuelling.py instead of re-writing the equations explicitly in constraint.py I get an inf error on the first iteration, would you be able to look at this as we had problems before about the constraint value

grmtrkngtn

I have a few questions about how this will handle a beam-fusion reactor, as this introduces some additional complications.

First, the thermal fuel mix is no longer necessarily 50/50, and the beam introduces both an additional fuel source and an additional fusion sink. I have left some comments in PlasmaFuelling on how the beam-target contribution could be accounted for in the species balances.

My understanding is that the existing PROCESS composition logic starts from electron density, calculates the total fuel-ion density, and then derives the individual D, T and He3 densities from the prescribed fuel fractions. If that is still the case, the composition routine may overwrite or constrain the same species densities that these new particle-balance constraints are intended to solve.

Have you also updated the density closure so that the absolute D and T densities can vary independently, with electron density and the fuel fractions then derived from charge neutrality? Otherwise, I am not sure that the D and T balances can act as independent solution constraints.

grmtrkngtn · 2026-07-22T08:32:24Z

+        )
+
+    @staticmethod
+    def calculate_deuterium_burnup_fraction(


There's an issue for beam target fusion tokamaks here, like VNS.

In a D beam case, the reaction is D_beam + T_thermal. Therefore thermal tritium is consumed but not thermal deuterium.

In this case, we cannot use the total DT rate, as it would overestimate D_thermal consumption, and therefore the required D fuelling.

Could we separate thermal DT from beam-target DT and create a beam-target sink according to the beam isotope fraction?

thermal_d_dt_consumption = (fusrat_dt_thermal + f_beam_tritium * fusrat_dt_beam)

Something like the above.

This will require a bit more work and most likely other constraints for the rate of fast D and T thermalisation. A bit like what will be added for the fast alphas. At the moment the total amount of D and T assumptions is just measured

grmtrkngtn · 2026-07-22T08:39:54Z

+        return 2 * fusrat_total / molflow_plasma_fuelling_vv_injected
+
+    @staticmethod
+    def calculate_tritium_burnup_fraction(


Also here we need to make a modification for the beam-fusion/VNS case.

For a D beam, both thermal–thermal DT and beam-target DT consume thermal tritium. For a T beam, the beam-target reaction consumes thermal deuterium instead.

Have added the fuelling components for the tritium beam

grmtrkngtn · 2026-07-22T08:43:54Z

+                * eta_plasma_fuelling
+                * molflow_plasma_fuelling_vv_injected
+            )
+            + fusrat_plasma_dhe3


Should this be a subtraction and not an addition? We're fuelling with He3 and consuming with fusion.

grmtrkngtn · 2026-07-28T14:31:57Z

    )

    # Deuterium and tritium ion densities
    nd_plasma_deuterium = nd_plasma_fuel_ions_vol_avg * f_deuterium_plasma


Should we add these to output? Might be useful/interesting.

We do have nd_plasma_fuel_ions_vol_avg which is the total fuel ion mix but not the individual species. May put a PR up to add this in a extra

…ng to be used to solve the fuelling and burnup equations

…riton, needed to know the indvidual fuel production and removal rates

…ude the equations for main fuel species flow and thermal alpha flow

…PlasmaFuelling model

…Physics class run workflow

…euterium, helium-3, and alpha particles

… for the burnup and fuelling rate now being calculated implicitly

…to previous functions and tests

…ass. Add details for some fusion reactions also

…ay in summary file

…resents

…itium source and loss rates. Implement these new methods in the constraints

Co-authored-by: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com>

…uterium source and loss rates

…g class; update constraints and flow rate calculations accordingly.

… update related constraints and documentation accordingly.

… thermal alpha particle source and loss rates; update related calculations in constraints and plotting functions.

…a class

…uelling and update related constraints and plotting functions

…e consistent notation for time derivatives.

… equilibrium solution models

…ak input file

chris-ashe added Physics Relating to the physics models Variable rename Input/Output Files Issues related to the input and output data files Documentation Improvements or additions to documentation Plotting labels Jun 29, 2026

chris-ashe assigned jonmaddock Jun 29, 2026

chris-ashe requested review from grmtrkngtn, j-a-foster, jonmaddock and timothy-nunn June 29, 2026 12:35

chris-ashe mentioned this pull request Jun 29, 2026

Particle balance constraint #4139

Closed

6 tasks

chris-ashe marked this pull request as ready for review June 29, 2026 12:36

chris-ashe requested a review from a team as a code owner June 29, 2026 12:36

j-a-foster requested changes Jun 29, 2026

View reviewed changes

jonmaddock requested changes Jun 30, 2026

View reviewed changes

chris-ashe force-pushed the particle_balance_constraint_fixed branch from 5c17efc to 8becc55 Compare July 21, 2026 08:48

grmtrkngtn requested changes Jul 22, 2026

View reviewed changes

chris-ashe force-pushed the particle_balance_constraint_fixed branch 4 times, most recently from 742f08b to 0d319b3 Compare July 23, 2026 12:54

grmtrkngtn reviewed Jul 28, 2026

View reviewed changes

Comment thread process/models/physics/physics.py Outdated

grmtrkngtn reviewed Jul 28, 2026

View reviewed changes

Comment thread process/models/physics/physics.py Outdated

grmtrkngtn reviewed Jul 28, 2026

View reviewed changes

Comment thread process/models/physics/physics.py Outdated

chris-ashe added 4 commits July 28, 2026 16:12

Add new input and iteration variables for plasma fuelling and recycli…

614521f

…ng to be used to solve the fuelling and burnup equations

Add Avogadro's number constant with reference documentation

dfc07cd

Add fusion rate density calculations for D-3He, D-D Helion, and D-D T…

f41e3f5

…riton, needed to know the indvidual fuel production and removal rates

Add PlasmaFuelling model class and integrate into Physics class. Incl…

3f93e55

…ude the equations for main fuel species flow and thermal alpha flow

chris-ashe and others added 25 commits July 28, 2026 16:12

Add fusion rate calculations for D-T and D-D reactions in Physics model

610561c

Add burnup fraction calculations for fuel, tritium, and deuterium in …

2452ee1

…PlasmaFuelling model

Implement run method for PlasmFuelling class and implement into main …

ad2cb4c

…Physics class run workflow

Add PlasmaFuelling model integration into Models class

30ebbc2

Add particle balance consistency equations constraints for tritium, d…

9d33c65

…euterium, helium-3, and alpha particles

Refactor phyaux method and related calculations in Physics to account…

8287fd0

… for the burnup and fuelling rate now being calculated implicitly

Remove now obsolete variables for phyaux and implement new variables …

028c87a

…to previous functions and tests

Add output method for fuelling info and refactor output in Physics cl…

448945b

…ass. Add details for some fusion reactions also

Calculate plasma fuel burnup fraction in Stellarator model

41e90fe

Add plasma fuelling documentation and update navigation in mkdocs

4b86c39

Add fuelling flow contour plots and update fuelling information displ…

336d354

…ay in summary file

Add new constraints to large_tokamak input and fix some runtime bugs

2d68d53

Add more descriptive notes to each function about what the output rep…

3437704

…resents

Refactor plasma fuelling calculations to use dedicated methods for tr…

f32b15a

…itium source and loss rates. Implement these new methods in the constraints

Update documentation/source/physics-models/plasma_fuelling.md

1fc3dfe

Co-authored-by: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com>

Update documentation/source/physics-models/plasma_fuelling.md

8261f1f

Co-authored-by: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com>

Quick fixes to correct for new thermal naming of alphas

fa58068

Refactor plasma fuelling calculations to use dedicated methods for de…

12f763d

…uterium source and loss rates

Implement helium-3 source and loss rate calculations in PlasmaFuellin…

a855acb

…g class; update constraints and flow rate calculations accordingly.

Refactor alpha particle flow calculations to use thermal terminology;…

44817ac

… update related constraints and documentation accordingly.

Refactor plasma fuelling calculations to include separate methods for…

8c54c4c

… thermal alpha particle source and loss rates; update related calculations in constraints and plotting functions.

Add beam deuterium and tritium injection rate variables to PhysicsDat…

7cbd2eb

…a class

Add beam deuterium and tritium injection rate calculations to PlasmaF…

02194e9

…uelling and update related constraints and plotting functions

Fix particle balance equations in Plasma Fuelling documentation to us…

634266c

…e consistent notation for time derivatives.

Update plasma fuelling documentation to clarify constraints as plasma…

2f8db51

… equilibrium solution models

chris-ashe force-pushed the particle_balance_constraint_fixed branch from 0d319b3 to 2f8db51 Compare July 28, 2026 15:12

chris-ashe added 3 commits July 28, 2026 16:14

Fix particle recycling and fuelling efficiency indices in large tokam…

a404f78

…ak input file

Remove duplicate output for reaction rates

1531d2b

Fix neutron production total calculation in Physics model

166d73a

Uh oh!

Conversation

chris-ashe commented Jun 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Checklist

Uh oh!

j-a-foster 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!

jonmaddock left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

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!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

chris-ashe commented Jul 1, 2026

Uh oh!

grmtrkngtn left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

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!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

chris-ashe commented Jun 29, 2026 •

edited

Loading