Particle balance constraint updated

[chris-ashe]

This PR add 5 new constraint equation with regard to plasma fuelling and particles loss constituting a way to enforce a plasma particle number balance for all species (except impurities)

To do so we now have our plasma fuelling into the vessel as an input parameter whose composition can be changed. The efficiency of this fuelling can be changed also to be indicative of the physical efficiency of the injection.

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.

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:

---

Solution changes

Before

After

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]

Why are these /sec and not /s?

[j-a-foster]

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

[jonmaddock]

• 1. 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?
• 2. Have you removed the original molflow_plasma_fuelling_required as an output? I might have missed that
• 3. 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?
• 4. 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?
• 5. How does this compare to the existing burnup calculation? Has it been removed?
• 6. 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?
• 7. I'm not sure about the fuelling composition constraint and total fuelling rate. Would individual species rates reduce the dimensionality?
• 8. 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).

The created docs were excellent. Thanks for the fixed-up commits, this was much easier to review.

[jonmaddock]

Not multiplied by 2 here.

[jonmaddock]

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

[chris-ashe]

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

[jonmaddock]

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?

[jonmaddock]

Snippet of constraints and opt params required to enable this please.

[chris-ashe]

@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]

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.

[chris-ashe]

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]

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.

[chris-ashe]

Have added the fuelling components for the tritium beam

[grmtrkngtn]

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

[grmtrkngtn]

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

[chris-ashe]

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

[jonmaddock]

Questions not yet answered, changes made unclear due to force-push.

[codecov-commenter]

Codecov Report

❌ Patch coverage is 44.19476% with 149 lines in your changes missing coverage. Please review.
✅ Project coverage is 49.39%. Comparing base (0be23b9) to head (d540f7d).

Files with missing lines	Patch %	Lines
process/core/io/plot/summary.py	6.74%	83 Missing ⚠️
process/models/physics/fuelling.py	54.21%	38 Missing ⚠️
process/models/physics/physics.py	4.76%	20 Missing ⚠️
process/core/solver/constraints.py	83.33%	4 Missing ⚠️
process/models/physics/fusion_reactions.py	0.00%	3 Missing ⚠️
process/models/stellarator/stellarator.py	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #4395      +/-   ##
==========================================
- Coverage   49.42%   49.39%   -0.03%     
==========================================
  Files         150      151       +1     
  Lines       30042    30278     +236     
==========================================
+ Hits        14847    14957     +110     
- Misses      15195    15321     +126     

☔ 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.

[chris-ashe]

I've numbered my questions now to aid the discussion.

1. Not done
2. Not sure if this has been answered in the docs (your 4th bullet) as I can't see the changes, but can you answer this here?
3. I think the "before and after" plots in the description partly answer this. Can you make clear what has changed in the "after" case (i.e. which constraints and optimisation parameters were added)? Given that the solution doesn't really change much, can you describe what is changing in the description (i.e. the fuelling rate, as in your 2nd bullet). Can you explain when these changes would have an effect on the solution (e.g. a machine with more asymmetric fuelling), and include this is the description?
4. This has been answered in your 3rd bullet, but has the old burnup calculation been removed? That's a fair increase! Is it right/agrees with known values?
5. Not done
6. Maybe in docs, but should be made clear here. Not done

I don't understand your 1st bullet: it's not clear to me what the system of equations now is. I think, as we discussed, it might be worth making a separate eval file that includes these modifications, so we have a "general" problem definition as we currently have and a separate "individual species" problem.

2. Yes, molflow_plasma_fuelling_required is now gone and replaced with molflow_plasma_fuelling_vv_injected, which is now our actual control optimisation parameter. So we give it as an input instead of letting it be calculated.
3. In the large tokamak optimisation case the final optimal solution has not changed. For the optimisation cases the number of equality constraints active has increased by 4. This is for the Tritium, Deuterium, Thermal Alpha, and fuelling compisition consistency. This is the same for the eval regression files. For the optimisation :

n_equality_constraints = 7

* Tritium particle balance
icc = 93

* Deuterium particle balance
icc = 94

* Alpha particle balance
icc = 96

* Fuelling composition consistency
icc = 97


----------------


* Particle recycling fraction
ixc = 178
f_plasma_particles_lcfs_recycled = 0.9

* Plasma fuelling efficiecy
ixc = 179
eta_plasma_fuelling = 0.7

* Injected VV fuelling rate
ixc = 180
molflow_plasma_fuelling_vv_injected = 5e21
boundl(180) = 1e20

* Deuterium fuelling fraction
ixc = 181
f_molflow_plasma_fuelling_deuterium = 0.5
boundl(181) = 0.4

* Tritium fuelling fraction
ixc = 182
f_molflow_plasma_fuelling_tritium = 0.5
boundl(182) = 0.4

We have the fuelling mixture, fuelling rate, fuelling efficiency and recycling coefficient as iteration variables.
For the eval case:

n_equality_constraints = 6

* Tritium particle balance
icc = 93
* Deuterium particle balance
icc = 94
* Alpha particle balance
icc = 96
* Fuelling composition consistency
icc = 97

-----------


* Particles recycled at LCFS
ixc = 178
f_plasma_particles_lcfs_recycled = 0.89

* Injected VV fuelling rate
ixc = 180
molflow_plasma_fuelling_vv_injected = 4.9992093e21

* Deuterium fuelling fraction
ixc = 181
f_molflow_plasma_fuelling_deuterium = 0.50092526

* Tritium fuelling fraction
ixc = 182
f_molflow_plasma_fuelling_tritium = 0.49907474

We have the same except we have the fuelling efficiency input so that our new number of iterations vars matches our increase increase in constraints (+4).

4. Yes the old burnup calc is now gone as we now calculate it implicitly based on its actual definition. Our old function used a guess via using the thermal alpha density, total ion density and the fusion rate with tauratio to make a loss term. This is not truly correct as you need to know the rate of creation vs the rate of injection. Our constraints now ensure that we know the true required rate of injection for all species. Therefore we can now calculate the equation implicitly as we already have the fusion rate for each species. The paper by Jackson "AN ANALYTIC EXPRESSION FOR THE TRITIUM
   BURNUP FRACTION IN BURNING-PLASMA DEVICES" tries to use the fuelling efficiency and recycling fraction but makes the use of parametrised profile values (that only work for parabolic) and other assumptions. In its verbose form for equation (1) he creates a reduced form of the equation he originally states which we actually do fully:
   .

In terms of the value for the burnup this can be heavily swayed by the recycling coefficient $R$ which really determines the system throughput. We could try and match the same value of $R$ to get the same burnup as before and it will tell us what our assumed recycling fractions would need to have been for our previous solutions

I will update the PR text to be more concise and give a problem setup for the input files @jonmaddock

[jonmaddock]

Ok, thanks @chris-ashe . So,

1. Done
2. Done
3. This is much clearer, thanks, but can you be more explicit? Include the full list of constraints as well as optimisation parameters for the optimisation and "eval" (solution) cases? Can you comment on why the fuelling efficiency is allowed to vary in the optimisation case but not in the solution case, and what it means (pellet to puff ratio, I think you said)?
4. Thanks for including the before and after plots, but I can't see any change in solution. Perhaps you could comment on what the difference is (however small), and under what circumstances these additional constraints and optimisation parameters would give a more different solution to the status quo?
5. (Burnup): Ok, does this mean that existing fuel ion equilibrium-only constraint doesn't work? Do we always need all of these new constraints on to make this new burnup calculation correct?
6. Is this done?
7. Done
8. Not done

Yours pedantically,

Jon

[chris-ashe]

Ok, thanks @chris-ashe . So,

1. Done
2. Done
3. This is much clearer, thanks, but can you be more explicit? Include the full list of constraints as well as optimisation parameters for the optimisation and "eval" (solution) cases? Can you comment on why the fuelling efficiency is allowed to vary in the optimisation case but not in the solution case, and what it means (pellet to puff ratio, I think you said)?
4. Thanks for including the before and after plots, but I can't see any change in solution. Perhaps you could comment on what the difference is (however small), and under what circumstances these additional constraints and optimisation parameters would give a more different solution to the status quo?
5. (Burnup): Ok, does this mean that existing fuel ion equilibrium-only constraint doesn't work? Do we always need all of these new constraints on to make this new burnup calculation correct?
6. Is this done?
7. Done
8. Not done

Yours pedantically,

Jon

3. I have ammended my previous comment to now show the constraints and total constraint number. The _eval case requires equal constraints and iterations variables, while in _optimisation case does not. The fuelling efficiency is allowed to change just to ease finding a solution. The problem is very open at the moment as their is no tritium plant or proper vacuum constraints applied to the fuelling inlet

4. For the _optimisation case their is no change to the solution for large_tokamak at the moment as the problem is open ended, meaning that the tritium plant and vacuum pumps arent putting a min/max on our fuelling inlet and outlet. When these are added we will probably see a change in fusion power.

5. Not sure what you mean by "mean that existing fuel ion equilibrium-only constraint doesn't work". For the burnup, since we know have individual species balances we now also have individual fuel species burnups (see output) these culminate together to make the total fuel burnup fraction. You only need the constraints on for the fuel species present in the plasma as it will enforce in=out. In most cases the 3He will not be needed as everyone just runs D-T. All current equations look for an exact analytical equation for the burnup which dosent exist due to the burnup problem requiring more than 1 ODE to solve for the full plasma case.

6. Please inspect the bottom of fuelling.md which contains a problem set up guide.

7. For the recycling fraction ($R$) and the fuelling efficiency the former is a higher fidelity unknown and the latter is a machine configuration facsimile. The recycling fraction normally has to be done by higher fidelity plasma-neutral modeling codes that couple the FW, Divertor and plasma interaction, we cannot do this here. The user will have to assume what their value is for their solution. Treat this the same as making the alpha to energy confinement ratio being equal to 5, its a prior higher order assumption. The fuelling efficiency will be based on what fuelling technology is used. This again is a 3D higher order number. What this means is we can now say "we require a fuelling efficiency > than ... etc"""

[jonmaddock]

3. Great. But

The fuelling efficiency is allowed to change just to ease finding a solution
Isn't a great reason. What does the fuelling efficiency actually mean? Is this something that can actually be changed at will, or is it an unknown? Otherwise this is a similar argument to having H-factor as an optimisation parameter.

4. Great, thanks. Done.
5. What I mean is does the existing fuel ion equilibrium constraint work by itself as it used to, using the original burnup fraction calculation? We may not necessarily want to enforce individual species equilibrium when trying to find an initial solution, as this will quite likely make it harder to converge. It also removes uncertain parameters (e.g. recycling fraction, fuelling efficiency) that a user may not know accurately/want to specify at an early stage of a design.
6. Great, done.
7. (Should the recycling fraction and fuelling efficiency be optimisation parameters?). I think this means that R and fuelling efficiency should be fixed, not used as parameters to solve these additional equations. There is a value of R, we just don't know it (form HiFi codes): that doesn't mean we get to pick its value to make our equations work and optimise.

[chris-ashe]

3. Great. But

The fuelling efficiency is allowed to change just to ease finding a solution
Isn't a great reason. What does the fuelling efficiency actually mean? Is this something that can actually be changed at will, or is it an unknown? Otherwise this is a similar argument to having H-factor as an optimisation parameter.

4. Great, thanks. Done.
5. What I mean is does the existing fuel ion equilibrium constraint work by itself as it used to, using the original burnup fraction calculation? We may not necessarily want to enforce individual species equilibrium when trying to find an initial solution, as this will quite likely make it harder to converge. It also removes uncertain parameters (e.g. recycling fraction, fuelling efficiency) that a user may not know accurately/want to specify at an early stage of a design.
6. Great, done.
7. (Should the recycling fraction and fuelling efficiency be optimisation parameters?). I think this means that R and fuelling efficiency should be fixed, not used as parameters to solve these additional equations. There is a value of R, we just don't know it (form HiFi codes): that doesn't mean we get to pick its value to make our equations work and optimise.

3. The fuelling efficiency is necessary component of the equation as in reality no plasma fuelling system is 100% perfect. Gas puffing is around 10% efficient and pellets are around up to 70%. For modelling this means we will always have a waste amount of fuel that will not penetrate the plasma and will need to be recycled in the divertor. This gives us realistic demand requirements for the vacuum pumps when they are properly integrated. In terms of it being an optimisation, allowing it to change can be seen as a facsimile for changing your fuelling method mix on the machine. Fuelling efficiency itself is too complex for us to model. The second benefit of having it is that the model can say we need a minimum efficiency to find solutions for an equilibrium, as a lower efficiency means more demand on the vacuum pumps which have their own limits. So we can say something like: "We need a fuelling efficiency of at least X in order to work with pump type Y. We currently have turbomolecular and cryopumps in our vacuum files which all have different pumping speeds for different species.

4. Which constraint do you mean exactly?

5. Again, this is similar to our alpha confinement ratio of 5. Realistically the relative fuelling efficiency can be changed in the control room via mass flow into the different systems, eg gas puffing and in pellets. The relative efficiency can also depend on the fuelling location control such as where on the plasma surface to puff or where in the core to aim the pellet. In terms of the recycling coefficient, most machines work in high recycling ($R = 0.9 - 0.99$) where the FW and divertor tiles are already full saturated so an escaped particle strikes the FW or divertor and is neutralised and returns to the plasma rather than being absorbed or pumped away. Advanced pumping with cryopumps can get $R$ down to around 0.8.