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 :

* 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:

* 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