From 2351cc53a6358ee5d5d8b74400a349499801fa19 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:11:37 +0100
Subject: [PATCH 01/39] Add new input and iteration variables for plasma
 fuelling and recycling to be used to solve the fuelling and burnup equations

---
 process/core/input.py                       | 16 ++++++
 process/core/solver/iteration_variables.py  | 26 ++++++++++
 process/data_structure/physics_variables.py | 57 +++++++++++++++++++++
 3 files changed, 99 insertions(+)

diff --git a/process/core/input.py b/process/core/input.py
index fb282b26ec..0008d256bc 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1182,6 +1182,22 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "f_len_sol_power_decay_inboard_outboard": InputVariable(
         "physics", float, range=(0.01, 2.0)
     ),
+    "f_plasma_particles_lcfs_recycled": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "eta_plasma_fuelling": InputVariable("physics", float, range=(0.0, 1.0)),
+    "molflow_plasma_fuelling_vv_injected": InputVariable(
+        "physics", float, range=(1e18, 1e24)
+    ),
+    "f_molflow_plasma_fuelling_deuterium": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "f_molflow_plasma_fuelling_tritium": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "f_molflow_plasma_fuelling_helium3": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
 }
 
 
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index d1df7c6feb..65c25b7f4b 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -239,6 +239,32 @@ class IterationVariable:
     175: IterationVariable("kappa", "physics", 0.00, 10.00),
     176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30),
     177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0),
+    178: IterationVariable("f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0),
+    179: IterationVariable("eta_plasma_fuelling", "physics", 0.01, 1.0),
+    180: IterationVariable(
+        "molflow_plasma_fuelling_vv_injected",
+        "physics",
+        1e18,
+        1e22,
+    ),
+    181: IterationVariable(
+        "f_molflow_plasma_fuelling_deuterium",
+        "physics",
+        0.0,
+        1.0,
+    ),
+    182: IterationVariable(
+        "f_molflow_plasma_fuelling_tritium",
+        "physics",
+        0.0,
+        1.0,
+    ),
+    183: IterationVariable(
+        "f_molflow_plasma_fuelling_helium3",
+        "physics",
+        0.0,
+        1.0,
+    ),
 }
 
 
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 53d7bf343a..6263a6a275 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -646,6 +646,15 @@ class PhysicsData:
     burnup: float = 0.0
     """fractional plasma burnup"""
 
+    f_plasma_fuel_burnup: float = 0.0
+    """Total fuel burnup fraction in plasma"""
+
+    f_plasma_tritium_burnup: float = 0.0
+    """Tritium burnup fraction in plasma"""
+
+    f_plasma_deuterium_burnup: float = 0.0
+    """Deuterium burnup fraction in plasma"""
+
     burnup_in: float = 0.0
     """fractional plasma burnup user input"""
 
@@ -803,6 +812,27 @@ class PhysicsData:
     fusrat_total: float = 0.0
     """fusion reaction rate, from beams and plasma (reactions/sec)"""
 
+    fusrat_plasma_dt: float = 0.0
+    """ D-T fusion reaction rate in plasma, (reactions/sec)"""
+
+    fusrat_dt_total: float = 0.0
+    """ Total D-T fusion reaction rate from beams and plasma, (reactions/sec)"""
+
+    fusrat_plasma_dd_helion: float = 0.0
+    """D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)"""
+
+    fusrat_plasma_dd_triton: float = 0.0
+    """D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)"""
+
+    fusrat_plasma_dd_total: float = 0.0
+    """Total D-D fusion reaction rate in plasma, (reactions/sec)"""
+
+    fusrat_plasma_dhe3: float = 0.0
+    """D-3He fusion reaction rate in plasma, (reactions/sec)"""
+
+    fusrat_neutron_production_total: float = 0.0
+    """Total neutron production rate from plasma and beams (neutrons/sec)"""
+
     fusrat_plasma_dt_profile: list[float] = field(default_factory=list)
     """Profile of D-T fusion reaction rate in plasma, (reactions/sec)"""
 
@@ -1326,6 +1356,33 @@ class PhysicsData:
     molflow_plasma_fuelling_required: float = 0.0
     """plasma fuelling rate (nucleus-pairs/s)"""
 
+    f_plasma_particles_lcfs_recycled: float = 0.9
+    """Fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)"""
+
+    eta_plasma_fuelling: float = 0.7
+    """Fuelling efficiency (η), fraction of fuel particles injected that become confined in the plasma"""
+
+    molflow_plasma_fuelling_vv_injected: float = 5e21
+    """Plasma fuelling rate into the vacuum vessel [particles/s]"""
+
+    molflow_plasma_fuelling_vv_injected_moles: float = 0.0
+    """Plasma fuelling rate into the vacuum vessel [moles/s]"""
+
+    molflow_plasma_fuelling_loss: float = 0.0
+    """Plasma fuelling rate that doesn't make it to plasma [particles/s]"""
+
+    molflow_plasma_fuelling_loss_moles: float = 0.0
+    """Plasma fuelling rate that doesn't make it to plasma [moles/s]"""
+
+    f_molflow_plasma_fuelling_deuterium: float = 0.5
+    """Fraction of plasma fuelling that is deuterium"""
+
+    f_molflow_plasma_fuelling_tritium: float = 0.5
+    """Fraction of plasma fuelling that is tritium"""
+
+    f_molflow_plasma_fuelling_helium3: float = 0.0
+    """Fraction of plasma fuelling that is helium-3"""
+
     tauratio: float = 1.0
     """tauratio /1.0/ : ratio of He and pellet particle confinement times"""
 

From 617d9316ab14588edc69fc9ae9ec9fbf4695cb33 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:12:20 +0100
Subject: [PATCH 02/39] Add Avogadro's number constant with reference
 documentation

---
 process/core/constants.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/process/core/constants.py b/process/core/constants.py
index 8583458d35..d6b8699de9 100644
--- a/process/core/constants.py
+++ b/process/core/constants.py
@@ -9,6 +9,12 @@
 MFILE = 13
 """Machine-optimised output file unit"""
 
+AVOGADRO_NUMBER = 6.02214076e23
+"""Avogadro's number [1/mol]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?na|search_for=avogadro
+"""
+
 ELECTRON_CHARGE = 1.602176634e-19
 """Electron / elementary charge [C]
 Reference: National Institute of Standards and Technology (NIST)

From 92b9a5caaea2d38f353d3815c1378b9c47c2a432 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:15:02 +0100
Subject: [PATCH 03/39] Add fusion rate density calculations for D-3He, D-D
 Helion, and D-D Triton, needed to know the indvidual fuel production and
 removal rates

---
 process/models/physics/fusion_reactions.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py
index 53d44266ff..5cceb0fcc0 100644
--- a/process/models/physics/fusion_reactions.py
+++ b/process/models/physics/fusion_reactions.py
@@ -347,6 +347,10 @@ def dhe3_reaction(self):
         alpha_rate_density = fusion_rate_density
         proton_rate_density = fusion_rate_density  # Proton production rate [m^3/second]
 
+        self.data.physics.fusrat_plasma_dhe3 = (
+            fusion_rate_density * self.data.physics.vol_plasma
+        )
+
         # Update the cumulative D-3He power density
         self.dhe3_power_density = fusion_power_density
 
@@ -451,6 +455,10 @@ def dd_helion_reaction(self):
         alpha_rate_density = 0.0
         proton_rate_density = 0.0
 
+        self.data.physics.fusrat_plasma_dd_helion = (
+            fusion_rate_density * self.data.physics.vol_plasma
+        )
+
         # Update the cumulative D-D power density
         self.dd_power_density += fusion_power_density
 
@@ -551,6 +559,10 @@ def dd_triton_reaction(self):
         # Proton production rate [particles/m³/s]
         proton_rate_density = fusion_rate_density
 
+        self.data.physics.fusrat_plasma_dd_triton = (
+            fusion_rate_density * self.data.physics.vol_plasma
+        )
+
         # Update the cumulative D-D power density
         self.dd_power_density += fusion_power_density
 

From e2abdc76ee0841a65ff64b4e1102081e9be45374 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:19:23 +0100
Subject: [PATCH 04/39] Add PlasmaFuelling model class and integrate into
 Physics class. Include the equations for main fuel species flow and thermal
 alpha flow

---
 process/models/physics/fuelling.py | 236 +++++++++++++++++++++++++++++
 process/models/physics/physics.py  |   3 +
 2 files changed, 239 insertions(+)
 create mode 100644 process/models/physics/fuelling.py

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
new file mode 100644
index 0000000000..e7d02acf88
--- /dev/null
+++ b/process/models/physics/fuelling.py
@@ -0,0 +1,236 @@
+import logging
+
+from process.core import constants
+from process.core.model import Model
+
+logger = logging.getLogger(__name__)
+
+
+class PlasmaFuelling(Model):
+    """Class to hold plasma fuelling calculations for plasma processing."""
+
+    def __init__(self):
+
+        self.outfile = constants.NOUT
+        self.mfile = constants.MFILE
+
+    def output(self):
+        """This model doesn't output to the output file, but it does generate contour
+        plots of plasma fuel flow rates vs recycling and fuelling efficiency.
+        """
+
+    @staticmethod
+    def calculate_plasma_tritium_flow_rate(
+        f_molflow_plasma_fuelling_tritium: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        fusrat_dt_total: float,
+        fusrat_plasma_dd_triton: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_tritium: float,
+    ) -> float:
+        """Calculate the tritium flow rate in the plasma exhaust.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_tritium : float
+            Fraction of tritium in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        fusrat_plasma_dd_triton : float
+            Tritium production rate from DD fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+        vol_plasma : float
+            Plasma volume (m^3).
+        f_plasma_fuel_tritium : float
+            Fraction of tritium in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Tritium flow rate in the plasma exhaust (particles/s).
+
+        """
+        return (
+            (
+                f_molflow_plasma_fuelling_tritium
+                * eta_plasma_fuelling
+                * molflow_plasma_fuelling_vv_injected
+            )
+            - fusrat_dt_total
+            + fusrat_plasma_dd_triton
+            - (
+                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium)
+                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
+            )
+        )
+
+    @staticmethod
+    def calculate_plasma_deuterium_flow_rate(
+        f_molflow_plasma_fuelling_deuterium: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        fusrat_dt_total: float,
+        fusrat_plasma_dhe3: float,
+        fusrat_plasma_dd_total: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_deuterium: float,
+    ) -> float:
+        """Calculate the deuterium flow rate in the plasma exhaust.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_deuterium : float
+            Fraction of deuterium in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+        fusrat_plasma_dd_total : float
+            Total deuterium consumption rate from DD fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+        vol_plasma : float
+            Plasma volume (m^3).
+        f_plasma_fuel_deuterium : float
+            Fraction of deuterium in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Deuterium flow rate in the plasma exhaust (particles/s).
+
+
+        """
+        return (
+            (
+                f_molflow_plasma_fuelling_deuterium
+                * eta_plasma_fuelling
+                * molflow_plasma_fuelling_vv_injected
+            )
+            - fusrat_dt_total
+            - 2 * fusrat_plasma_dd_total
+            - fusrat_plasma_dhe3
+            - (
+                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium)
+                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
+            )
+        )
+
+    @staticmethod
+    def calculate_plasma_helium3_flow_rate(
+        f_molflow_plasma_fuelling_helium3: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        fusrat_plasma_dhe3: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_helium3: float,
+    ) -> float:
+        """Calculate the helium-3 flow rate in the plasma exhaust.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_helium3 : float
+            Fraction of helium-3 in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+        vol_plasma : float
+            Plasma volume (m^3).
+        f_plasma_fuel_helium3 : float
+            Fraction of helium-3 in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Helium-3 flow rate in the plasma exhaust (particles/s).
+
+        """
+        return (
+            (
+                f_molflow_plasma_fuelling_helium3
+                * eta_plasma_fuelling
+                * molflow_plasma_fuelling_vv_injected
+            )
+            + fusrat_plasma_dhe3
+            - (
+                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_helium3)
+                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
+            )
+        )
+
+    @staticmethod
+    def calculate_plasma_alphas_flow_rate(
+        fusrat_dt_total: float,
+        fusrat_plasma_dhe3: float,
+        t_energy_confinement: float,
+        f_t_alpha_energy_confinement: float,
+        nd_plasma_alphas_vol_avg: float,
+        vol_plasma: float,
+    ) -> float:
+        """Calculate the alpha particle flow rate in the plasma exhaust.
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_t_alpha_energy_confinement : float
+            Ratio of alpha particle confinement time to energy confinement time (dimensionless).
+        nd_plasma_alphas_vol_avg : float
+            Volume-averaged density of alpha particles in the plasma (particles/m^3).
+        vol_plasma : float
+            Plasma volume (m^3).
+
+        Returns
+        -------
+        float
+            Alpha particle flow rate in the plasma exhaust (particles/s).
+
+        """
+        # Alpha particle balance
+
+        return (
+            fusrat_dt_total
+            + fusrat_plasma_dhe3
+            - (nd_plasma_alphas_vol_avg * vol_plasma)
+            / (t_energy_confinement * f_t_alpha_energy_confinement)
+        )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index b6e1ce43dc..bf591a9de6 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -39,6 +39,7 @@
     )
     from process.models.physics.density_limit import PlasmaDensityLimit
     from process.models.physics.exhaust import PlasmaExhaust
+    from process.models.physics.fuelling import PlasmaFuelling
     from process.models.physics.l_h_transition import PlasmaConfinementTransition
     from process.models.physics.plasma_current import (
         PlasmaCurrent,
@@ -199,6 +200,7 @@ def __init__(
         plasma_dia_current: PlasmaDiamagneticCurrent,
         plasma_geometry: PlasmaGeom,
         scrape_off_layer: ScrapeOffLayer,
+        plasma_fuelling: PlasmaFuelling,
     ):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
@@ -216,6 +218,7 @@ def __init__(
         self.dia_current = plasma_dia_current
         self.geometry = plasma_geometry
         self.scrape_off_layer = scrape_off_layer
+        self.fuelling = plasma_fuelling
 
     def output(self) -> None:
         """Output plasma physics information."""

From 6e3e8de39b26c2de88c3fa13481ce0c2222cce89 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:20:34 +0100
Subject: [PATCH 05/39] Add fusion rate calculations for D-T and D-D reactions
 in Physics model

---
 process/models/physics/physics.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index bf591a9de6..b2c78a6a77 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -683,6 +683,22 @@ def run(self):
             self.data.physics.fusden_total * self.data.physics.vol_plasma
         )
 
+        self.data.physics.fusrat_plasma_dt = (self.data.physics.p_plasma_dt_mw * 1e6) / (
+            constants.D_T_ENERGY
+        )
+        self.data.physics.fusrat_plasma_dd_total = (
+            self.data.physics.fusrat_plasma_dd_helion
+            + self.data.physics.fusrat_plasma_dd_triton
+        )
+
+        self.data.physics.fusrat_neutron_production_total = (
+            self.data.physics.fusrat_plasma_dd_helion + self.data.physics.fusrat_dt_total
+        )
+
+        self.data.physics.fusrat_dt_total = (
+            self.data.physics.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY)
+        )
+
         # Create some derived values and add beam contribution to fusion power
         (
             self.data.physics.pden_neutron_total_mw,

From 1e8cf30aa9ab03f3a526f6068436ec4189b1639f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:21:29 +0100
Subject: [PATCH 06/39] Add burnup fraction calculations for fuel, tritium, and
 deuterium in PlasmaFuelling model

---
 process/models/physics/fuelling.py | 93 ++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index e7d02acf88..ad4dc6000e 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -19,6 +19,99 @@ def output(self):
         plots of plasma fuel flow rates vs recycling and fuelling efficiency.
         """
 
+    @staticmethod
+    def calculate_fuel_burnup_fraction(
+        fusrat_total: float, molflow_plasma_fuelling_vv_injected: float
+    ) -> float:
+        """Calculate the fuel burnup fraction
+
+        Parameters
+        ----------
+        fusrat_total : float
+            Total fusion rate (particles/s).
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate into vacuum vessel (particles/s).
+
+        Returns
+        -------
+        Fuel burnup fraction (dimensionless).
+
+        Notes
+        -----
+        The fusion rate is multiplied by two to convert from nucleus pairs to particles,
+        as the fuelling rate is in particles/s.
+
+        """
+        return 2 * fusrat_total / molflow_plasma_fuelling_vv_injected
+
+    @staticmethod
+    def calculate_tritium_burnup_fraction(
+        fusrat_dt_total: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        f_molflow_plasma_fuelling_tritium: float,
+    ) -> float:
+        """Calculate the tritium burnup fraction
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate into vacuum vessel (particles/s).
+        f_molflow_plasma_fuelling_tritium : float
+            Fraction of tritium in the plasma fuelling.
+
+        Returns
+        -------
+        Tritium burnup fraction (dimensionless).
+
+        Notes
+        -----
+        The fusion rate is multiplied by two to convert from nucleus pairs to particles,
+        as the fuelling rate is in particles/s.
+
+        """
+        return fusrat_dt_total / (
+            molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_tritium
+        )
+
+    @staticmethod
+    def calculate_deuterium_burnup_fraction(
+        fusrat_dt_total: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        f_molflow_plasma_fuelling_deuterium: float,
+        fusrat_plasma_dd_total: float,
+        fusrat_plasma_dhe3: float,
+    ) -> float:
+        """Calculate the deuterium burnup fraction
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate into vacuum vessel (particles/s).
+        f_molflow_plasma_fuelling_deuterium : float
+            Fraction of deuterium in the plasma fuelling.
+        fusrat_plasma_dd_total : float
+            Total deuterium consumption rate from DD fusion (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+
+        Returns
+        -------
+        Deuterium burnup fraction (dimensionless).
+
+        Notes
+        -----
+        The fusion rate is multiplied by two to convert from nucleus pairs to particles,
+        as the fuelling rate is in particles/s.
+
+        """
+        return (fusrat_dt_total + 2 * fusrat_plasma_dd_total + fusrat_plasma_dhe3) / (
+            molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium
+        )
+
     @staticmethod
     def calculate_plasma_tritium_flow_rate(
         f_molflow_plasma_fuelling_tritium: float,

From 395623e17e249558879dd4d74e5a9bd7a3ce2450 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:23:06 +0100
Subject: [PATCH 07/39] Implement run method for PlasmFuelling class and
 implement into main Physics class run workflow

---
 process/models/physics/fuelling.py | 32 ++++++++++++++++++++++++++++++
 process/models/physics/physics.py  |  2 ++
 2 files changed, 34 insertions(+)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index ad4dc6000e..1e4409a93b 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -14,6 +14,38 @@ def __init__(self):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
 
+    def run(self):
+        self.data.physics.molflow_plasma_fuelling_vv_injected_moles = (
+            self.data.physics.molflow_plasma_fuelling_vv_injected
+            / constants.AVOGADRO_NUMBER
+        )
+
+        self.data.physics.molflow_plasma_fuelling_loss = (
+            self.data.physics.molflow_plasma_fuelling_vv_injected
+            * (1 - self.data.physics.eta_plasma_fuelling)
+        )
+        self.data.physics.molflow_plasma_fuelling_loss_moles = (
+            self.data.physics.molflow_plasma_fuelling_loss / constants.AVOGADRO_NUMBER
+        )
+
+        self.data.physics.f_plasma_fuel_burnup = self.calculate_fuel_burnup_fraction(
+            fusrat_total=self.data.physics.fusrat_total,
+            molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+        )
+        self.data.physics.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction(
+            fusrat_dt_total=self.data.physics.fusrat_dt_total,
+            molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+            f_molflow_plasma_fuelling_tritium=self.data.physics.f_molflow_plasma_fuelling_tritium,
+        )
+
+        self.data.physics.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction(
+            fusrat_plasma_dd_total=self.data.physics.fusrat_plasma_dd_total,
+            molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+            f_molflow_plasma_fuelling_deuterium=self.data.physics.f_molflow_plasma_fuelling_deuterium,
+            fusrat_dt_total=self.data.physics.fusrat_dt_total,
+            fusrat_plasma_dhe3=self.data.physics.fusrat_plasma_dhe3,
+        )
+
     def output(self):
         """This model doesn't output to the output file, but it does generate contour
         plots of plasma fuel flow rates vs recycling and fuelling efficiency.
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index b2c78a6a77..3bb48f7fbb 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -998,6 +998,8 @@ def run(self):
             self.data.physics.tauratio,
         )
 
+        self.fuelling.run()
+
         self.data.physics.ntau, self.data.physics.nTtau = (
             self.confinement.calculate_double_and_triple_product(
                 nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg,

From 2adfcfba6e4249e9a3d0dbe08b3dad33786c9542 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:24:34 +0100
Subject: [PATCH 08/39] Add PlasmaFuelling model integration into Models class

---
 process/main.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/process/main.py b/process/main.py
index 4f2c7ccfa8..4e282de5ff 100644
--- a/process/main.py
+++ b/process/main.py
@@ -80,6 +80,7 @@
 )
 from process.models.physics.density_limit import PlasmaDensityLimit
 from process.models.physics.exhaust import PlasmaExhaust
+from process.models.physics.fuelling import PlasmaFuelling
 from process.models.physics.impurity_radiation import (
     initialise_imprad,
 )
@@ -693,6 +694,7 @@ def __init__(self, data: DataStructure):
         self.plasma_density_limit = PlasmaDensityLimit()
         self.plasma_exhaust = PlasmaExhaust()
         self.sauter_bootstrap_current = SauterBootstrapCurrent()
+        self.plasma_fuelling = PlasmaFuelling()
         self.plasma_bootstrap_current = PlasmaBootstrapCurrent(
             plasma_profile=self.plasma_profile,
             sauter_bootstrap=self.sauter_bootstrap_current,
@@ -718,6 +720,7 @@ def __init__(self, data: DataStructure):
             plasma_dia_current=self.plasma_dia_current,
             plasma_geometry=self.plasma_geom,
             scrape_off_layer=self.scrape_off_layer,
+            plasma_fuelling=self.plasma_fuelling,
         )
         self.physics_detailed = DetailedPhysics(
             plasma_profile=self.plasma_profile,
@@ -824,6 +827,7 @@ def models(self) -> tuple[Model, ...]:
             self.plasma_fields,
             self.sauter_bootstrap_current,
             self.plasma_transition,
+            self.self.plasma_fuelling,
             self.physics_detailed,
             self.electron_cyclotron,
             self.lower_hybrid,

From ac144dad00d100c5b2f1376722829cca48698ff0 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:26:39 +0100
Subject: [PATCH 09/39] Add particle balance consistency equations constraints
 for tritium, deuterium, helium-3, and alpha particles

---
 process/core/solver/constraints.py | 147 +++++++++++++++++++++++++++++
 process/data_structure/numerics.py |   5 +
 2 files changed, 152 insertions(+)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 72c11014fd..5b68c9194a 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -1955,6 +1955,153 @@ def constraint_equation_92(constraint_registration, data):
     )
 
 
+@ConstraintManager.register_constraint(93, "particles/s", "=")
+def constraint_equation_93(constraint_registration, data):
+    """
+    Tritium particle balance consistency equation.
+
+    This constraint ensures that the tritium particle balance is maintained in the
+    plasma.
+
+    The numerator represents the total tritium sources, including fuelling and fusion
+    reactions, while the denominator represents the total tritium sinks, which are
+    related to the plasma volume, fuel ion density, and confinement time.
+
+    The equation ensures that the rate of tritium production equals the rate of tritium
+    loss, maintaining a steady-state condition for tritium in the plasma.
+    """
+    numerator = (
+        (
+            data.physics.eta_plasma_fuelling
+            * data.physics.molflow_plasma_fuelling_vv_injected
+            * data.physics.f_molflow_plasma_fuelling_tritium
+        )
+        + data.physics.fusrat_plasma_dd_triton
+        - data.physics.fusrat_dt_total
+    )
+    denominator = (
+        data.physics.nd_plasma_fuel_ions_vol_avg
+        * data.physics.vol_plasma
+        * data.physics.f_plasma_fuel_tritium
+    ) / (
+        data.physics.t_energy_confinement
+        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    )
+
+    return eq(numerator, denominator, constraint_registration)
+
+
+@ConstraintManager.register_constraint(94, "particles/s", "=")
+def constraint_equation_94(constraint_registration, data):
+    """
+    Deuterium particle balance consistency equation.
+
+    This constraint ensures that the deuterium particle balance is maintained in the
+    plasma.
+
+    The numerator represents the total deuterium sources, including fuelling and fusion
+    reactions, while the denominator represents the total deuterium sinks, which are
+    related to the plasma volume, fuel ion density, and confinement time.
+
+    The equation ensures that the rate of deuterium production equals the rate of deuterium
+    loss, maintaining a steady-state condition for deuterium in the plasma.
+    """
+    numerator = (
+        (
+            data.physics.eta_plasma_fuelling
+            * data.physics.molflow_plasma_fuelling_vv_injected
+            * data.physics.f_molflow_plasma_fuelling_deuterium
+        )
+        - data.physics.fusrat_dt_total
+        - data.physics.fusrat_plasma_dhe3
+        - 2.0 * data.physics.fusrat_plasma_dd_total
+    )
+    denominator = (
+        data.physics.nd_plasma_fuel_ions_vol_avg
+        * data.physics.vol_plasma
+        * data.physics.f_plasma_fuel_deuterium
+    ) / (
+        data.physics.t_energy_confinement
+        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    )
+
+    return eq(numerator, denominator, constraint_registration)
+
+
+@ConstraintManager.register_constraint(95, "particles/s", "=")
+def constraint_equation_95(constraint_registration, data):
+    """
+    Helium-3 particle balance consistency equation.
+
+    This constraint ensures that the helium-3 particle balance is maintained in the
+    plasma.
+
+    The numerator represents the total helium-3 sources, including fuelling and fusion
+    reactions, while the denominator represents the total helium-3 sinks, which are
+    related to the plasma volume, fuel ion density, and confinement time.
+
+    The equation ensures that the rate of helium-3 production equals the rate of helium-3
+    loss, maintaining a steady-state condition for helium-3 in the plasma.
+    """
+    numerator = (
+        data.physics.eta_plasma_fuelling
+        * data.physics.molflow_plasma_fuelling_vv_injected
+        * data.physics.f_molflow_plasma_fuelling_helium3
+    ) + data.physics.fusrat_plasma_dhe3
+    denominator = (
+        data.physics.nd_plasma_fuel_ions_vol_avg
+        * data.physics.vol_plasma
+        * data.physics.f_plasma_fuel_helium3
+    ) / (
+        data.physics.t_energy_confinement
+        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    )
+
+    return eq(numerator, denominator, constraint_registration)
+
+
+@ConstraintManager.register_constraint(96, "particles/s", "=")
+def constraint_equation_96(constraint_registration, data):
+    """
+    Alpha particle balance consistency equation.
+
+    This constraint ensures that the alpha particle balance is maintained in the
+    plasma.
+
+    The numerator represents the total alpha particle sources,
+    while the denominator represents the total alpha particle sinks, which are
+    related to the plasma volume, fuel ion density, and confinement time.
+
+    The equation ensures that the rate of alpha particle production equals the rate of
+    alpha particle loss, maintaining a steady-state condition for alpha particles in
+    the plasma.
+    """
+    # Alpha particle balance
+    numerator = data.physics.fusrat_dt_total + data.physics.fusrat_plasma_dhe3
+    denominator = (data.physics.nd_plasma_alphas_vol_avg * data.physics.vol_plasma) / (
+        data.physics.t_energy_confinement * data.physics.f_alpha_energy_confinement
+    )
+
+    return eq(numerator, denominator, constraint_registration)
+
+
+@ConstraintManager.register_constraint(97, "", "=")
+def constraint_equation_97(constraint_registration, data):
+    """Equation for checking the fuelling composition is consistent.
+
+    f_molflow_plasma_fuelling_deuterium: fraction of deuterium ions
+    f_molflow_plasma_fuelling_tritium: fraction of tritium ions
+    f_molflow_plasma_fuelling_helium3: fraction of helium-3 ions
+    """
+    return eq(
+        data.physics.f_molflow_plasma_fuelling_deuterium
+        + data.physics.f_molflow_plasma_fuelling_tritium
+        + data.physics.f_molflow_plasma_fuelling_helium3,
+        1.0,
+        constraint_registration,
+    )
+
+
 def constraint_eqns(m: int, ieqn: int, data: DataStructure):
     """Evaluates the constraints given the current state of PROCESS.
 
diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py
index 43b0b78314..1c306b767b 100644
--- a/process/data_structure/numerics.py
+++ b/process/data_structure/numerics.py
@@ -292,6 +292,11 @@ class NumericsData:
             "CS achievable stress load cycles lower limit           ",
             "ECRH ignitability                ",  # Stellarator constraint
             "Fuel composition consistency     ",
+            "Tritium particle balance consistency        ",
+            "Deuterium particle balance consistency     ",
+            "Helium-3 particle balance consistency     ",
+            "Alpha particle balance consistency     ",
+            "Fuelling composition consistency     ",
         ]
     )
     """Labels describing constraint equations (corresponding itvs)<UL>

From 1529860c55aee571287a84f6ee03b718c1e3a890 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:36:12 +0100
Subject: [PATCH 10/39] Refactor phyaux method and related calculations in
 Physics to account for the burnup and fuelling rate now being calculated
 implicitly

---
 process/models/physics/physics.py         | 91 ++---------------------
 process/models/stellarator/stellarator.py | 14 ----
 tests/unit/models/physics/test_physics.py | 91 -----------------------
 3 files changed, 5 insertions(+), 191 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 3bb48f7fbb..ce3329ffd9 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -974,28 +974,13 @@ def run(self):
             0.5e0 * self.data.physics.ind_plasma * self.data.physics.plasma_current**2
         )
 
-        # Calculate auxiliary physics related information
-        sbar = 1.0e0
         (
-            self.data.physics.burnup,
-            self.data.physics.figmer,
-            self.data.physics.fusrat,
-            self.data.physics.molflow_plasma_fuelling_required,
-            self.data.physics.rndfuel,
             self.data.physics.t_alpha_confinement,
             self.data.physics.f_t_alpha_energy_confinement,
         ) = self.phyaux(
-            self.data.physics.aspect,
-            self.data.physics.nd_plasma_fuel_ions_vol_avg,
-            self.data.physics.fusden_total,
-            self.data.physics.fusden_alpha_total,
-            self.data.physics.plasma_current,
-            sbar,
-            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
-            self.data.physics.t_energy_confinement,
-            self.data.physics.vol_plasma,
-            self.data.physics.burnup_in,
-            self.data.physics.tauratio,
+            fusden_alpha_total=self.data.physics.fusden_alpha_total,
+            nd_plasma_alphas_vol_avg=self.data.physics.nd_plasma_alphas_thermal_vol_avg,
+            t_energy_confinement=self.data.physics.t_energy_confinement,
         )
 
         self.fuelling.run()
@@ -1504,66 +1489,32 @@ def plasma_composition(self):
     @staticmethod
     @nb.njit(cache=True)
     def phyaux(
-        aspect: float,
-        nd_plasma_fuel_ions_vol_avg: float,
-        fusden_total: float,
         fusden_alpha_total: float,
-        plasma_current: float,
-        sbar: float,
         nd_plasma_alphas_thermal_vol_avg: float,
         t_energy_confinement: float,
-        vol_plasma: float,
-        burnup_in: float,
-        tauratio: float,
-    ) -> tuple[float, float, float, float, float, float, float, float]:
+    ) -> tuple[float, float]:
         """Auxiliary physics quantities
 
         Parameters
         ----------
-        aspect : float
-            Plasma aspect ratio.
-        nd_plasma_fuel_ions_vol_avg : float
-            Fuel ion density (/m3).
-        fusden_total : float
-            Fusion reaction rate from plasma and beams (/m3/s).
         fusden_alpha_total : float
             Alpha particle production rate (/m3/s).
-        plasma_current : float
-            Plasma current (A).
-        sbar : float
-            Exponent for aspect ratio (normally 1).
         nd_plasma_alphas_thermal_vol_avg : float
             Alpha ash density (/m3).
         t_energy_confinement : float
             Global energy confinement time (s).
-        vol_plasma : float
-            Plasma volume (m3).
-        burnup_in: float
-            fractional plasma burnup user input
-        tauratio: float
-            ratio of He and pellet particle confinement times
+
 
         Returns
         -------
         tuple
             A tuple containing:
-            - burnup (float): Fractional plasma burnup.
-            - figmer (float): Physics figure of merit.
-            - fusrat (float): Number of fusion reactions per second.
-            - molflow_plasma_fuelling_required (float): Fuelling rate for D-T
-              (nucleus-pairs/sec).
-            - rndfuel (float): Fuel burnup rate (reactions/s).
             - t_alpha_confinement (float): Alpha particle confinement time (s).
             - f_t_alpha_energy_confinement (float): Fraction of alpha energy confinement.
             This subroutine calculates extra physics related items needed by other
             parts of the code.
 
         """
-        figmer = 1e-6 * plasma_current * aspect**sbar
-
-        # Fusion reactions per second
-        fusrat = fusden_total * vol_plasma
-
         # Alpha particle confinement time (s)
         # Number of alphas / alpha production rate
         # only likely if DD is only active fusion reaction
@@ -1573,41 +1524,9 @@ def phyaux(
             else nd_plasma_alphas_thermal_vol_avg / fusden_alpha_total
         )
 
-        # Fractional burnup
-        # (Consider detailed model in: G. L. Jackson, V. S. Chan, R. D. Stambaugh,
-        # Fusion Science and Technology, vol.64, no.1, July 2013, pp.8-12)
-        # The ratio of ash to fuel particle confinement times is given by
-        # tauratio
-        # Possible logic...
-        # burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3;
-        # fuel ion-pairs burned/m3 = alpha particles/m3 (for both D-T and
-        # D-He3 reactions)
-        # initial fuel ion-pairs/m3 = burnt fuel ion-pairs/m3 + unburnt fuel-ion
-        # pairs/m3
-        # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3
-        if burnup_in <= 1.0e-9:
-            burnup = (
-                nd_plasma_alphas_thermal_vol_avg
-                / (nd_plasma_alphas_thermal_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg)
-                / tauratio
-            )
-        else:
-            burnup = burnup_in
-
-        # Fuel burnup rate (reactions/second) (previously Amps)
-        rndfuel = fusrat
-
-        # Required fuelling rate (fuel ion pairs/second) (previously Amps)
-        molflow_plasma_fuelling_required = rndfuel / burnup
-
         f_t_alpha_energy_confinement = t_alpha_confinement / t_energy_confinement
 
         return (
-            burnup,
-            figmer,
-            fusrat,
-            molflow_plasma_fuelling_required,
-            rndfuel,
             t_alpha_confinement,
             f_t_alpha_energy_confinement,
         )
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index e75fc106c8..6804db193f 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2376,27 +2376,13 @@ def st_phys(self, output):
         # Calculate auxiliary physics related information
         # for the rest of the code
 
-        sbar = 1.0e0
         (
-            self.data.physics.burnup,
-            self.data.physics.figmer,
-            _fusrat,
-            self.data.physics.molflow_plasma_fuelling_required,
-            self.data.physics.rndfuel,
             self.data.physics.t_alpha_confinement,
             self.data.physics.f_t_alpha_energy_confinement,
         ) = self.physics.phyaux(
-            self.data.physics.aspect,
-            self.data.physics.nd_plasma_fuel_ions_vol_avg,
-            self.data.physics.fusden_total,
             self.data.physics.fusden_alpha_total,
-            self.data.physics.plasma_current,
-            sbar,
             self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             self.data.physics.t_energy_confinement,
-            self.data.physics.vol_plasma,
-            self.data.physics.burnup_in,
-            self.data.physics.tauratio,
         )
 
         # Calculate the neoclassical sanity check with PROCESS parameters
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 9e3aa6f47b..cc578d7cfe 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1837,46 +1837,12 @@ def test_vscalc(voltsecondreqparam):
 
 
 class PhyauxParam(NamedTuple):
-    tauratio: Any = None
-
-    burnup_in: Any = None
-
-    aspect: Any = None
-
-    nd_plasma_electrons_vol_avg: Any = None
-
-    te: Any = None
-
-    nd_plasma_fuel_ions_vol_avg: Any = None
-
     nd_plasma_alphas_thermal_vol_avg: Any = None
 
-    fusden_total: Any = None
-
     fusden_alpha_total: Any = None
 
-    plasma_current: Any = None
-
-    sbar: Any = None
-
     t_energy_confinement: Any = None
 
-    vol_plasma: Any = None
-
-    expected_burnup: Any = None
-
-    expected_ntau: Any = None
-
-    expected_nTtau: Any = None
-
-    expected_figmer: Any = None
-
-    expected_fusrat: Any = None
-
-    expected_molflow_plasma_fuelling_required: Any = None
-
-    expected_rndfuel: Any = None
-
     expected_t_alpha_confinement: Any = None
 
 
@@ -1884,41 +1850,13 @@ class PhyauxParam(NamedTuple):
     "phyauxparam",
     [
         PhyauxParam(
-            tauratio=1,
-            burnup_in=0,
-            aspect=3,
-            nd_plasma_fuel_ions_vol_avg=5.858890125e19,
-            nd_plasma_alphas_thermal_vol_avg=7.5e18,
-            fusden_total=1.9852091609123786e17,
             fusden_alpha_total=1.973996644759543e17,
-            plasma_current=18398455.678867526,
-            sbar=1,
             t_energy_confinement=3.401323521525641,
-            vol_plasma=1888.1711539956691,
-            expected_burnup=0.20383508579699033,
-            expected_figmer=55.195367036602576,
-            expected_fusrat=3.7484146722826997e20,
-            expected_molflow_plasma_fuelling_required=1.838944781084418e21,
-            expected_rndfuel=3.7484146722826997e20,
             expected_t_alpha_confinement=37.993985551650177,
         ),
         PhyauxParam(
-            tauratio=1,
-            burnup_in=0,
-            aspect=3,
-            nd_plasma_fuel_ions_vol_avg=5.8576156204039725e19,
-            nd_plasma_alphas_thermal_vol_avg=7.5e18,
-            fusden_total=1.9843269653375773e17,
             fusden_alpha_total=1.9731194318497056e17,
-            plasma_current=18398455.678867526,
-            sbar=1,
             t_energy_confinement=3.402116961408892,
-            vol_plasma=1888.1711539956691,
-            expected_burnup=0.20387039462081086,
-            expected_figmer=55.195367036602576,
-            expected_fusrat=3.7467489360461772e20,
-            expected_molflow_plasma_fuelling_required=1.8378092331723546e21,
-            expected_rndfuel=3.7467489360461772e20,
             expected_t_alpha_confinement=38.010876984618747,
         ),
     ],
@@ -1937,44 +1875,15 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
     :type monkeypatch: _pytest.monkeypatch.monkeypatch
     """
 
-    monkeypatch.setattr(physics.data.physics, "tauratio", phyauxparam.tauratio)
-
-    monkeypatch.setattr(physics.data.physics, "burnup_in", phyauxparam.burnup_in)
-
     (
-        burnup,
-        figmer,
-        fusrat,
-        molflow_plasma_fuelling_required,
-        rndfuel,
         t_alpha_confinement,
         _,
     ) = physics.phyaux(
-        aspect=phyauxparam.aspect,
-        nd_plasma_fuel_ions_vol_avg=phyauxparam.nd_plasma_fuel_ions_vol_avg,
         nd_plasma_alphas_thermal_vol_avg=phyauxparam.nd_plasma_alphas_thermal_vol_avg,
-        fusden_total=phyauxparam.fusden_total,
         fusden_alpha_total=phyauxparam.fusden_alpha_total,
-        plasma_current=phyauxparam.plasma_current,
-        sbar=phyauxparam.sbar,
         t_energy_confinement=phyauxparam.t_energy_confinement,
-        vol_plasma=phyauxparam.vol_plasma,
-        burnup_in=phyauxparam.burnup_in,
-        tauratio=phyauxparam.tauratio,
     )
 
-    assert burnup == pytest.approx(phyauxparam.expected_burnup)
-
-    assert figmer == pytest.approx(phyauxparam.expected_figmer)
-
-    assert fusrat == pytest.approx(phyauxparam.expected_fusrat)
-
-    assert molflow_plasma_fuelling_required == pytest.approx(
-        phyauxparam.expected_molflow_plasma_fuelling_required
-    )
-
-    assert rndfuel == pytest.approx(phyauxparam.expected_rndfuel)
-
     assert t_alpha_confinement == pytest.approx(phyauxparam.expected_t_alpha_confinement)
 
 

From a88831ce2674cb8c5f7ce77d81f2867a75237bb7 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:42:13 +0100
Subject: [PATCH 11/39] Remove now obsolete variables for phyaux and implement
 new variables to previous functions and tests

---
 process/data_structure/physics_variables.py | 12 ------------
 process/models/costs/costs.py               |  2 +-
 process/models/vacuum.py                    |  8 ++++----
 tests/unit/models/test_costs_1990.py        | 12 +++++++-----
 tests/unit/models/test_vacuum.py            |  2 +-
 5 files changed, 13 insertions(+), 23 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 6263a6a275..a51ad2325e 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -643,9 +643,6 @@ class PhysicsData:
     e_plasma_magnetic_stored: float = 0.0
     """Plasma stored magnetic energy [J]"""
 
-    burnup: float = 0.0
-    """fractional plasma burnup"""
-
     f_plasma_fuel_burnup: float = 0.0
     """Total fuel burnup fraction in plasma"""
 
@@ -1353,9 +1350,6 @@ class PhysicsData:
     """Plasma safety factor at 95% flux surface (q₉₅) (`iteration variable 18`)
     """
 
-    molflow_plasma_fuelling_required: float = 0.0
-    """plasma fuelling rate (nucleus-pairs/s)"""
-
     f_plasma_particles_lcfs_recycled: float = 0.9
     """Fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)"""
 
@@ -1383,9 +1377,6 @@ class PhysicsData:
     f_molflow_plasma_fuelling_helium3: float = 0.0
     """Fraction of plasma fuelling that is helium-3"""
 
-    tauratio: float = 1.0
-    """tauratio /1.0/ : ratio of He and pellet particle confinement times"""
-
     q95_min: float = 0.0
     """Plasmalower limit for edge safety factor"""
 
@@ -1431,9 +1422,6 @@ class PhysicsData:
     f_nd_plasma_carbon_electron: float = 0.0
     """n_carbon / n_e"""
 
-    rndfuel: float = 0.0
-    """fuel burnup rate (reactions/second)"""
-
     f_nd_plasma_iron_argon_electron: float = 0.0
     """n_highZ / n_e"""
 
diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index 208f98cbe0..2de33737f6 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -2354,7 +2354,7 @@ def acc2272(self):
             #  New calculation: 2 nuclei * reactions/sec * kg/nucleus * g/kg * sec/day
             self.data.physics.wtgpd = (
                 2.0e0
-                * self.data.physics.rndfuel
+                * self.data.physics.fusrat_total
                 * self.data.physics.m_fuel_amu
                 * constants.UMASS
                 * 1000.0e0
diff --git a/process/models/vacuum.py b/process/models/vacuum.py
index b305bd93ea..ea1b929d27 100644
--- a/process/models/vacuum.py
+++ b/process/models/vacuum.py
@@ -49,7 +49,7 @@ def run(self, output: bool = False):
         # MDK Check this!!
         gasld = (
             2.0e0
-            * self.data.physics.molflow_plasma_fuelling_required
+            * self.data.physics.molflow_plasma_fuelling_vv_injected
             * self.data.physics.m_fuel_amu
             * constants.UMASS
         )
@@ -112,7 +112,7 @@ def vacuum_simple(self, output) -> float:
         # One ITER torus cryopump has a throughput of 50 Pa m3/s = 1.2155e+22 molecules/s
         # Issue #304
         n_iter_vacuum_pumps = (
-            self.data.physics.molflow_plasma_fuelling_required
+            self.data.physics.molflow_plasma_fuelling_vv_injected
             / self.data.vacuum.molflow_vac_pumps
         )
 
@@ -161,8 +161,8 @@ def _vacuum_simple_output(self, n_iter_vacuum_pumps, npumpdown, npump):
         process_output.ovarre(
             self.outfile,
             "Plasma fuelling rate (nucleus-pairs/s)",
-            "(molflow_plasma_fuelling_required)",
-            self.data.physics.molflow_plasma_fuelling_required,
+            "(molflow_plasma_fuelling_vv_injected)",
+            self.data.physics.molflow_plasma_fuelling_vv_injected,
             "OP ",
         )
 
diff --git a/tests/unit/models/test_costs_1990.py b/tests/unit/models/test_costs_1990.py
index 6d442bc2c6..9095ef4cf5 100644
--- a/tests/unit/models/test_costs_1990.py
+++ b/tests/unit/models/test_costs_1990.py
@@ -151,7 +151,7 @@ def test_acc2272(monkeypatch, costs):
     :param monkeypatch: Mock fixture
     :type monkeypatch: object
     """
-    monkeypatch.setattr(costs.data.physics, "rndfuel", 7.158e20)
+    monkeypatch.setattr(costs.data.physics, "fusrat_total", 7.158e20)
     monkeypatch.setattr(costs.data.physics, "m_fuel_amu", 2.5)
     monkeypatch.setattr(costs.data.costs, "fkind", 1)
     monkeypatch.setattr(costs.data.costs, "c2271", 0)
@@ -4031,7 +4031,9 @@ class Acc2272Param(NamedTuple):
 
     edrive: Any = None
 
-    rndfuel: Any = None
+    wtgpd: Any = None
+
+    fusrat_total: Any = None
 
     m_fuel_amu: Any = None
 
@@ -4054,7 +4056,7 @@ class Acc2272Param(NamedTuple):
             ife=0,
             gain=0,
             edrive=5000000,
-            rndfuel=7.0799717510383796e20,
+            fusrat_total=7.0799717510383796e20,
             m_fuel_amu=2.5,
             c227=0,
             c22=0,
@@ -4068,7 +4070,7 @@ class Acc2272Param(NamedTuple):
             ife=0,
             gain=0,
             edrive=5000000,
-            rndfuel=7.0777619721108953e20,
+            fusrat_total=7.0777619721108953e20,
             m_fuel_amu=2.5,
             c227=284.96904049038437,
             c22=3474.7391916096453,
@@ -4103,7 +4105,7 @@ def test_acc2272_rut(acc2272param, monkeypatch, costs):
 
     monkeypatch.setattr(costs.data.ife, "edrive", acc2272param.edrive)
 
-    monkeypatch.setattr(costs.data.physics, "rndfuel", acc2272param.rndfuel)
+    monkeypatch.setattr(costs.data.physics, "fusrat_total", acc2272param.fusrat_total)
 
     monkeypatch.setattr(costs.data.physics, "m_fuel_amu", acc2272param.m_fuel_amu)
 
diff --git a/tests/unit/models/test_vacuum.py b/tests/unit/models/test_vacuum.py
index f1913200a2..db1bb461b0 100644
--- a/tests/unit/models/test_vacuum.py
+++ b/tests/unit/models/test_vacuum.py
@@ -40,7 +40,7 @@ def test_simple_model(monkeypatch, vacuum):
         """
         monkeypatch.setattr(
             vacuum.data.physics,
-            "molflow_plasma_fuelling_required",
+            "molflow_plasma_fuelling_vv_injected",
             7.5745668997694112e22,
         )
         monkeypatch.setattr(vacuum.data.physics, "a_plasma_surface", 1500.3146527709359)

From 23d849be81e09922ed87b1ae7af45f16c64ba3ee Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:46:38 +0100
Subject: [PATCH 12/39] Add output method for fuelling info and refactor output
 in Physics class. Add details for some fusion reactions also

---
 process/models/physics/fuelling.py | 129 ++++++++++++++++++++++++++++-
 process/models/physics/physics.py  | 122 ++++++++++++++++-----------
 2 files changed, 197 insertions(+), 54 deletions(-)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 1e4409a93b..b9e70799b8 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -1,9 +1,10 @@
-import logging
-
 from process.core import constants
+from process.core import process_output as po
 from process.core.model import Model
-
-logger = logging.getLogger(__name__)
+from process.data_structure import (
+    numerics,
+    reinke_variables,
+)
 
 
 class PlasmaFuelling(Model):
@@ -359,3 +360,123 @@ def calculate_plasma_alphas_flow_rate(
             - (nd_plasma_alphas_vol_avg * vol_plasma)
             / (t_energy_confinement * f_t_alpha_energy_confinement)
         )
+
+    def output_fuelling_info(self):
+        """Output fuelling information to mfile."""
+        po.oheadr(self.outfile, "Plasma Fuelling")
+        po.ovarre(
+            self.outfile,
+            "Fuelling rate (nucleus-pairs/s)",
+            "(molflow_plasma_fuelling_vv_injected)",
+            self.data.physics.molflow_plasma_fuelling_vv_injected,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fuelling rate (moles/s)",
+            "(molflow_plasma_fuelling_vv_injected_moles)",
+            self.data.physics.molflow_plasma_fuelling_vv_injected_moles,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fuelling loss (nucleus-pairs/s)",
+            "(molflow_plasma_fuelling_loss)",
+            self.data.physics.molflow_plasma_fuelling_loss,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fuelling loss (moles/s)",
+            "(molflow_plasma_fuelling_loss_moles)",
+            self.data.physics.molflow_plasma_fuelling_loss_moles,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Fraction of plasma fuelling that is deuterium",
+            "(f_molflow_plasma_fuelling_deuterium)",
+            self.data.physics.f_molflow_plasma_fuelling_deuterium,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fraction of plasma fuelling that is tritium",
+            "(f_molflow_plasma_fuelling_tritium)",
+            self.data.physics.f_molflow_plasma_fuelling_tritium,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fraction of plasma fuelling that is helium-3",
+            "(f_molflow_plasma_fuelling_helium3)",
+            self.data.physics.f_molflow_plasma_fuelling_helium3,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Fuelling efficiency",
+            "(eta_plasma_fuelling)",
+            self.data.physics.eta_plasma_fuelling,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fraction of plasma particles recycled at the LCFS",
+            "(f_plasma_particles_lcfs_recycled)",
+            self.data.physics.f_plasma_particles_lcfs_recycled,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fuel burn-up rate (reactions/s)",
+            "(fusrat_total)",
+            self.data.physics.fusrat_total,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+        po.ovarrf(
+            self.outfile,
+            "Total fuel burn-up fraction",
+            "(f_plasma_fuel_burnup)",
+            self.data.physics.f_plasma_fuel_burnup,
+            "OP ",
+        )
+        po.ovarrf(
+            self.outfile,
+            "Tritium burn-up fraction",
+            "(f_plasma_tritium_burnup)",
+            self.data.physics.f_plasma_tritium_burnup,
+            "OP ",
+        )
+        po.ovarrf(
+            self.outfile,
+            "Deuterium burn-up fraction",
+            "(f_plasma_deuterium_burnup)",
+            self.data.physics.f_plasma_deuterium_burnup,
+            "OP ",
+        )
+
+        if 78 in numerics.icc:
+            po.osubhd(self.outfile, "Reinke Criterion :")
+            po.ovarin(
+                self.outfile,
+                "index of impurity to be iterated for divertor detachment",
+                "(impvardiv)",
+                reinke_variables.impvardiv,
+            )
+            po.ovarre(
+                self.outfile,
+                "Minimum Impurity fraction from Reinke",
+                "(fzmin)",
+                reinke_variables.fzmin,
+                "OP ",
+            )
+            po.ovarre(
+                self.outfile,
+                "Actual Impurity fraction",
+                "(fzactual)",
+                reinke_variables.fzactual,
+            )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index ce3329ffd9..9b2fb36588 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -1805,6 +1805,77 @@ def outplas(self):
             self.data.physics.fusrat_total,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "D-T Fusion rate: total (reactions/sec)",
+            "(fusrat_dt_total)",
+            self.data.physics.fusrat_dt_total,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-T Fusion rate: plasma (reactions/sec)",
+            "(fusrat_plasma_dt)",
+            self.data.physics.fusrat_plasma_dt,
+            "OP ",
+        )
+
+        po.ovarre(
+            self.outfile,
+            "D-D -> 3He Fusion rate: plasma (reactions/sec)",
+            "(fusrat_plasma_dd_helion)",
+            self.data.physics.fusrat_plasma_dd_helion,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-D -> T Fusion rate: plasma (reactions/sec)",
+            "(fusrat_plasma_dd_triton)",
+            self.data.physics.fusrat_plasma_dd_triton,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-D Fusion rate: total (reactions/sec)",
+            "(fusrat_plasma_dd_total)",
+            self.data.physics.fusrat_plasma_dd_total,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-3He Fusion rate: total (reactions/sec)",
+            "(fusrat_plasma_dhe3)",
+            self.data.physics.fusrat_plasma_dhe3,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Neutron production rate: total (particles/sec)",
+            "(fusrat_neutron_production_total)",
+            self.data.physics.fusrat_neutron_production_total,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-3He Fusion rate: total (reactions/sec)",
+            "(fusrat_plasma_dhe3)",
+            self.data.physics.fusrat_plasma_dhe3,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Neutron production rate: total (particles/sec)",
+            "(fusrat_neutron_production_total)",
+            self.data.physics.fusrat_neutron_production_total,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "D-D Fusion rate: total (reactions/sec)",
+            "(fusrat_plasma_dd_total)",
+            self.data.physics.fusrat_plasma_dd_total,
+            "OP ",
+        )
         po.ovarre(
             self.outfile,
             "Fusion rate density: total (reactions/m³/sec)",
@@ -2648,56 +2719,7 @@ def outplas(self):
             self.plasma_bootstrap_current.output()
             self.dia_current.output()
 
-        po.oheadr(self.outfile, "Plasma Fuelling")
-        po.ovarre(
-            self.outfile,
-            "Ratio of He and pellet particle confinement times",
-            "(tauratio)",
-            self.data.physics.tauratio,
-        )
-        po.ovarre(
-            self.outfile,
-            "Fuelling rate (nucleus-pairs/s)",
-            "(molflow_plasma_fuelling_required)",
-            self.data.physics.molflow_plasma_fuelling_required,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Fuel burn-up rate (reactions/s)",
-            "(rndfuel)",
-            self.data.physics.rndfuel,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Burn-up fraction",
-            "(burnup)",
-            self.data.physics.burnup,
-            "OP ",
-        )
-
-        if 78 in self.data.numerics.icc:
-            po.osubhd(self.outfile, "Reinke Criterion :")
-            po.ovarre(
-                self.outfile,
-                "index of impurity to be iterated for divertor detachment",
-                "(impvardiv)",
-                self.data.reinke.impvardiv,
-            )
-            po.ovarre(
-                self.outfile,
-                "Minimum Impurity fraction from Reinke",
-                "(fzmin)",
-                self.data.reinke.fzmin,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Actual Impurity fraction",
-                "(fzactual)",
-                self.data.reinke.fzactual,
-            )
+        self.fuelling.output_fuelling_info()
 
     def output_temperature_density_profile_info(self) -> None:
         """Output information about plasma temperature and density profiles."""

From 395d81ec5723c682092926b5cc9cee8799f0d7d2 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:47:20 +0100
Subject: [PATCH 13/39] Calculate plasma fuel burnup fraction in Stellarator
 model

---
 process/models/stellarator/stellarator.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 6804db193f..05d004deb2 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2385,6 +2385,11 @@ def st_phys(self, output):
             self.data.physics.t_energy_confinement,
         )
 
+        self.data.physics.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction(
+            fusrat_total=self.data.physics.fusrat_total,
+            molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+        )
+
         # Calculate the neoclassical sanity check with PROCESS parameters
         (
             q_PROCESS,

From 281bfe68cdf784c37208a4a60923bf9867e23e12 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:49:28 +0100
Subject: [PATCH 14/39] Add plasma fuelling documentation and update navigation
 in mkdocs

---
 .../source/physics-models/plasma_fuelling.md  | 209 ++++++++++++++++++
 mkdocs.yml                                    |   1 +
 2 files changed, 210 insertions(+)
 create mode 100644 documentation/source/physics-models/plasma_fuelling.md

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
new file mode 100644
index 0000000000..b9a05389e3
--- /dev/null
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -0,0 +1,209 @@
+# Plasma Fuelling | `PlasmaFuelling()`
+
+## Particle balance
+
+The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$.
+
+$$
+\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+$$
+
+$$
+\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+$$
+
+$$
+\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+$$
+
+$$
+\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*}
+$$
+
+In a steady state equilibrium all 4 of these equations should balance, therefore:
+
+$$
+\frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha}}{dt} = 0
+$$
+
+Here $\eta_{\text{fuelling}}$ is the fuelling efficiecny which represents the method of injecting fuel into the plasma. Gas puffing on the low field side is probably around 0.01-0.1, supersonic gas is 0.1 and 0.2 and using pellets can get you close to unity with 0.5-0.9. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel, so $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ together presents the fraction of injected fuel that actually makes it into the plasma core to fuse. 
+
+
+The fuelling fractional compositions is given by $f$
+
+ - $N$ is the total amount of ions in the plasma.
+
+ - $\tau_{\text{fuel}}^*$ is the recycling corrected fuel particle confinement time given by:
+
+    $\tau_{\text{fuel}}^* = (\tau_p) /  (1-R)$
+
+The (effective) exhaust efficiency is is given by , $\eta_{\text{eff}} = 1- R$ 
+
+The factor $\frac{R}{1-R}$ is the mean number of recycling events back into the burning region experieced by a particle before it is pumped away.
+
+The definition of the recycling coefficient $R = 1- \frac{\Gamma_{\text{pumps}}}{\Gamma_{\text{out}}}$, where $\Gamma_{\text{pumps}}$ is the number of particles exhausted by the pumps per second and $\Gamma_{\text{out}}$ is the number of particles per second transported radially outwards across the separatrix.
+
+
+
+Where $\tau_p$ is the particle confinement time which we can assume is approximately equal to the energy confinement time ($\tau_p = \tau_E$). 
+
+!!! warning "Relation between $\tau_p$ and $\tau_E$"
+
+    In this model the "raw" particle confinement time ($\tau_p$) is set to always match the energy confinement time ($\tau_E$). The only variation of this in terms of the "effective" particle confinement time $\tau_{p}^*$ is via the recycling coefficient ($R$). If needed new coeffcients may be added to scale ($\tau_p$) before recyling corrections if needed.
+
+!!! note "Quantifying $R$"
+
+    The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1].
+
+
+
+
+### METIS Alpha Confinement
+
+$$
+\tau_{\alpha} = f_{\alpha}\tau_{\text{E}}\frac{R}{1-R}\tau_{\text{ne}}
+$$
+
+This is the model currently in METIS[^2] and is found in [^3]
+
+
+--------------
+
+### Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()`
+
+$$
+\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+$$
+
+---------------
+
+### Deuterium Flow Rate | `calculate_plasma_deuterium_flow_rate()`
+
+$$
+\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+$$
+
+---------------
+
+### Helium-3 Flow Rate | `calculate_plasma_helium3_flow_rate()`
+
+$$
+\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+$$
+
+---------------
+
+### Alpha Particle Flow Rate | `calculate_plasma_alphas_flow_rate()`
+
+$$
+\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*}
+$$
+
+-----------------
+
+## Fuel Burnup Fration
+
+In a steady state tokamak, the burnup fraction ($f_b$) is explicitly defined by the following rate equation:
+
+$$
+f_b = \frac{\Gamma_{\text{fusion}}}{\Gamma_{\text{fuel}}}
+$$
+
+where $\Gamma_{\text{fusion}}$ is the fusion rate in reactions per second $[\text{s}^{-1}]$, and $\Gamma_{\text{fuel}}$ is our fuelling rate into the vessel like is shown above also in $[\text{s}^{-1}]$.
+
+-----------------
+
+### Total Fuel Burnup Fraction | `calculate_fuel_burnup_fraction()`
+
+For the total burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_fuel_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+T}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}}
+$$
+
+Here the factor of 2 is included as each fusion reaction removes 2 particles but our fuelling rate looks at indivudal particles injected.
+
+-------------------
+
+### Tritium Burnup Fraction | `calculate_tritium_burnup_fraction()`
+
+For just the tritium burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_tritium_burnup}} = \frac{\left(\Gamma_{\text{D+T}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,T}}} 
+$$
+
+------------------
+
+### Deuterium Burnup Fraction | `calculate_deuterium_burnup_fraction()`
+
+For just the deuterium burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_deuterium_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,D}}}
+$$
+
+------------------
+
+## Key Constraints
+
+### Tritium Flow Consistency
+
+This constraint can be activated by stating `icc = 93` in the input file.
+
+This constraint ensures that the change in tritium particles as a function of time is zero. It ensures the output of `calculate_plasma_tritium_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma consistency model**
+
+-----------------
+
+### Deuterium Flow Consistency
+
+This constraint can be activated by stating `icc = 94` in the input file.
+
+This constraint ensures that the change in deuterium particles as a function of time is zero. It ensures the output of `calculate_plasma_deuterium_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma consistency model**
+
+----------------
+
+### Helium-3 Flow Consistency
+
+This constraint can be activated by stating `icc = 95` in the input file.
+
+This constraint ensures that the change in helium-3 particles as a function of time is zero. It ensures the output of `calculate_plasma_helium3_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma consistency model**
+
+----------------
+
+### Alpha Particle Flow Consistency
+
+This constraint can be activated by stating `icc = 96` in the input file.
+
+This constraint ensures that the change in alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_flow_rate(()` is zero
+
+**It is recommended to have this constraint on as it is a plasma consistency model**
+
+------------------
+
+### Fuelling Proportion Consistency
+
+This constraint can be activated by stating `icc = 97` in the input file.
+
+This ensures that all 3 injected fuelling fractions sum up to 1:
+
+$$
+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**
+
+-----------------
+
+
+[^1]: G. L. Jackson, V. S. Chan, and R. D. Stambaugh, “An Analytic Expression for the Tritium Burnup Fraction in Burning-Plasma Devices,” Fusion Science and Technology, vol. 64, no. 1, pp. 8–12, Jul. 2013, doi: https://doi.org/10.13182/fst13-a17042.
+
+[^2]: J. F. Artaud et al., “Metis: a fast integrated tokamak modelling tool for scenario design,” Nuclear Fusion, vol. 58, no. 10, pp. 105001–105001, Aug. 2018, doi: https://doi.org/10.1088/1741-4326/aad5b1.
+
+[^3]: D. Reiter, H. Kever, G. H. Wolf, M. Baelmans, R. Behrisch, and R. Schneider, “Helium removal from tokamaks,” Plasma Physics and Controlled Fusion, vol. 33, no. 13, pp. 1579–1600, Nov. 1991, doi: https://doi.org/10.1088/0741-3335/33/13/008.
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index 6cfa762788..cbdaf4def1 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -53,6 +53,7 @@ nav:
                   - Overview: physics-models/plasma_beta/plasma_beta.md
                   - Fast Alpha: physics-models/plasma_beta/plasma_alpha_beta_contribution.md
               - Density Limit: physics-models/plasma_density.md
+              - Fuelling: physics-models/plasma_fuelling.md
               - Composition & Impurities: physics-models/plasma_composition.md
               - Radiation: physics-models/plasma_radiation.md
               - Plasma Current:

From 1dedb527d1c3751b5f3a261df4ff3df13aa28c7b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 11:58:35 +0100
Subject: [PATCH 15/39] Add fuelling flow contour plots and update fuelling
 information display in summary file

---
 process/core/io/plot/summary.py | 399 +++++++++++++++++++++++++++++---
 1 file changed, 364 insertions(+), 35 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index bd34580d43..f903cd97b0 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -68,6 +68,7 @@
 )
 from process.models.physics.density_limit import DensityLimitModel
 from process.models.physics.exhaust import calculate_brunner_divertor_power_splits
+from process.models.physics.fuelling import PlasmaFuelling
 from process.models.physics.impurity_radiation import read_impurity_file
 from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel
 from process.models.physics.physics import (
@@ -2942,7 +2943,7 @@ def plot_main_plasma_information(
         f"|  D: {mfile.get('f_plasma_fuel_deuterium', scan=scan):.2f}  |  T: {mfile.get('f_plasma_fuel_tritium', scan=scan):.2f}  |  3He: {mfile.get('f_plasma_fuel_helium3', scan=scan):.2f}  |\n\n"
         f"Fusion Power, $P_{{\\text{{fus}}}}:$ {mfile.get('p_fusion_total_mw', scan=scan):,.2f} MW\n"
         f"D-T Power, $P_{{\\text{{fus,DT}}}}:$ {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n"
-        f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n"
+        f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n"
         f"D-3He Power, $P_{{\\text{{fus,D3He}}}}:$ {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW\n"
         f"Alpha Power, $P_{{\\alpha}}:$ {mfile.get('p_alpha_total_mw', scan=scan):,.2f} MW"
     )
@@ -2966,9 +2967,9 @@ def plot_main_plasma_information(
         f"   - Average mass of all plasma ions: {mfile.get('m_ions_total_amu', scan=scan):.3f} amu\n"
         f"Fuel mass: {mfile.get('m_plasma_fuel_ions', scan=scan) * 1000:.4f} g\n"
         f"   - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n"
-        f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n"
-        f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n"
-        f"Burn-up fraction: {mfile.get('burnup', scan=scan):.4f} \n"
+        f"Fueling rate: {mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.3e} nucleus-pairs/s\n"
+        f"Fuel burn-up rate: {mfile.get('fusrat_total', scan=scan):.3e} reactions/s \n"
+        f"Burn-up fraction: {mfile.get('f_plasma_fuel_burnup', scan=scan):.4f} \n"
     )
 
     axis.text(
@@ -12113,8 +12114,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     # Add plasma volume, areas and shaping information
     textstr_general = (
         f"Total fusion rate: {mfile.get('fusrat_total', scan=scan):.4e} reactions/s\n"
-        f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m3/s\n"
-        f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m3/s\n"
+        f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m$^3$/s\n"
+        f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m$^3$/s"
     )
 
     axis.text(
@@ -12136,9 +12137,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     # ============================================================================
 
     textstr_dt = (
-        f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n"
-        f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.2f} MW                     \n"
-        f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.2f} MW\n"
+        f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.4f} MW\n"
+        f"Total fusion rate: {mfile.get('fusrat_dt_total', scan=scan):.4e} reactions/s\n"
+        f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.4f} MW                              \n"
+        f"Plasma fusion rate: {mfile.get('fusrat_plasma_dt', scan=scan):.4e} reactions/s\n"
+        f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.4f} MW"
     )
 
     axis.text(
@@ -12158,7 +12161,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     )
 
     axis.text(
-        0.24,
+        0.285,
         0.8,
         "$\\text{D - T}$",
         fontsize=20,
@@ -12169,13 +12172,16 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     # =================================================
 
     textstr_dd = (
-        f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n"
-        f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f}                      \n"
+        f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n"
+        f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f}                                  \n\n"
+        f"D+D -> T fusion rate: {mfile.get('fusrat_plasma_dd_triton', scan=scan):.4e} reactions/s                       \n"
+        f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s                     \n"
+        f"Total D-D fusion rate: {mfile.get('fusrat_plasma_dd_total', scan=scan):.4e} reactions/s"
     )
 
     axis.text(
         0.05,
-        0.65,
+        0.625,
         textstr_dd,
         fontsize=9,
         verticalalignment="bottom",
@@ -12190,8 +12196,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     )
 
     axis.text(
-        0.22,
-        0.685,
+        0.31,
+        0.69,
         "$\\text{D - D}$",
         fontsize=20,
         verticalalignment="top",
@@ -12200,8 +12206,10 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
 
     # =================================================
 
-    textstr_dhe3 = f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW                                 \n\n"
-
+    textstr_dhe3 = (
+        f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.4f} MW                                         \n"
+        f"D+3He fusion rate: {mfile.get('fusrat_plasma_dhe3', scan=scan):.4e} reactions/s                               \n"
+    )
     axis.text(
         0.05,
         0.55,
@@ -12219,8 +12227,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     )
 
     axis.text(
-        0.21,
-        0.59,
+        0.285,
+        0.56,
         "$\\text{D - 3He}$",
         fontsize=20,
         verticalalignment="top",
@@ -12230,15 +12238,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     # =================================================
 
     textstr_alpha = (
-        f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.2f} MW\n"
-        f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.2f} MW\n"
-        f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.2f} MW\n\n"
-        f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m3/sec\n"
-        f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m3/sec\n\n"
-        f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n"
-        f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n"
-        f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n"
-        f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n"
+        f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.4f} MW\n"
+        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"
+        f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m$^3$\n"
+        f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m$^3$\n\n"
+        f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m$^3$\n"
+        f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m$^3$"
     )
 
     axis.text(
@@ -12269,11 +12277,12 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     # =================================================
 
     textstr_neutron = (
-        f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.2f} MW\n"
-        f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.2f} MW\n"
-        f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.2f} MW\n\n"
-        f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m3\n"
-        f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m3\n"
+        f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.4f} MW\n"
+        f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.4f} MW\n"
+        f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.4f} MW\n\n"
+        f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m$^3$\n"
+        f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m$^3$\n\n"
+        f"Neutron production rate: {mfile.get('fusrat_neutron_production_total', scan=scan):.4e} particles/s"
     )
 
     axis.text(
@@ -12293,8 +12302,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
     )
 
     axis.text(
-        0.25,
-        0.2,
+        0.26,
+        0.21,
         "$n$",
         fontsize=20,
         verticalalignment="top",
@@ -16474,6 +16483,311 @@ def plot_pf_dimensions(
     axis.set_aspect("equal", adjustable="box")
 
 
+def plot_tritium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
+    """Plot contour of tritium flow rate vs recycling and fuelling rate."""
+    recycling_range = np.linspace(0.01, 0.99, 20)
+    fuelling_range = np.linspace(0.01, 1.0, 20)
+    tritium_flow = np.zeros((len(recycling_range), len(fuelling_range)))
+
+    for i, recycling in enumerate(recycling_range):
+        for j, fuelling in enumerate(fuelling_range):
+            tritium_flow[i, j] = PlasmaFuelling.calculate_plasma_tritium_flow_rate(
+                f_molflow_plasma_fuelling_tritium=mfile.get(
+                    "f_molflow_plasma_fuelling_tritium", scan=scan
+                ),
+                eta_plasma_fuelling=fuelling,
+                molflow_plasma_fuelling_vv_injected=mfile.get(
+                    "molflow_plasma_fuelling_vv_injected", scan=scan
+                ),
+                fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan),
+                fusrat_plasma_dd_triton=mfile.get("fusrat_plasma_dd_triton", scan=scan),
+                t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
+                f_plasma_particles_lcfs_recycled=recycling,
+                nd_plasma_fuel_ions_vol_avg=mfile.get(
+                    "nd_plasma_fuel_ions_vol_avg", scan=scan
+                ),
+                vol_plasma=mfile.get("vol_plasma", scan=scan),
+                f_plasma_fuel_tritium=mfile.get("f_plasma_fuel_tritium", scan=scan),
+            )
+
+    contour = axis.contourf(
+        fuelling_range,
+        recycling_range,
+        tritium_flow,
+        levels=15,
+        cmap="seismic",
+        norm=plt.matplotlib.colors.CenteredNorm(vcenter=0),
+    )
+
+    axis.contour(
+        fuelling_range,
+        recycling_range,
+        tritium_flow,
+        levels=[0],
+        colors="black",
+        linewidths=2,
+    )
+
+    # Plot star for mfile values
+    recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan)
+    fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan)
+    axis.plot(
+        fuelling_mfile,
+        recycling_mfile,
+        marker="*",
+        markersize=15,
+        color="yellow",
+        markeredgecolor="black",
+        markeredgewidth=1.5,
+    )
+
+    axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)")
+    axis.set_ylabel("Recycling Fraction [$R$]")
+    axis.set_title("Plasma Tritium Flow Rate (particles/s)", pad=20)
+    axis.minorticks_on()
+    axis.grid(True, which="major", linestyle="-", alpha=0.7)
+    axis.grid(True, which="minor", linestyle=":", alpha=0.4)
+    cbar = plt.colorbar(contour, ax=axis, label="Tritium Flow Rate")
+    cbar.ax.axhline(y=0, color="black", linewidth=2)
+
+
+def plot_deuterium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
+    """Plot contour of deuterium flow rate vs recycling and fuelling rate."""
+    recycling_range = np.linspace(0.01, 0.99, 20)
+    fuelling_range = np.linspace(0.01, 1.0, 20)
+    deuterium_flow = np.zeros((len(recycling_range), len(fuelling_range)))
+
+    for i, recycling in enumerate(recycling_range):
+        for j, fuelling in enumerate(fuelling_range):
+            deuterium_flow[i, j] = PlasmaFuelling.calculate_plasma_deuterium_flow_rate(
+                f_molflow_plasma_fuelling_deuterium=mfile.get(
+                    "f_molflow_plasma_fuelling_deuterium", scan=scan
+                ),
+                eta_plasma_fuelling=fuelling,
+                molflow_plasma_fuelling_vv_injected=mfile.get(
+                    "molflow_plasma_fuelling_vv_injected", scan=scan
+                ),
+                fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan),
+                fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
+                fusrat_plasma_dd_total=mfile.get("fusrat_plasma_dd_total", scan=scan),
+                t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
+                f_plasma_particles_lcfs_recycled=recycling,
+                nd_plasma_fuel_ions_vol_avg=mfile.get(
+                    "nd_plasma_fuel_ions_vol_avg", scan=scan
+                ),
+                vol_plasma=mfile.get("vol_plasma", scan=scan),
+                f_plasma_fuel_deuterium=mfile.get("f_plasma_fuel_deuterium", scan=scan),
+            )
+
+    contour = axis.contourf(
+        fuelling_range,
+        recycling_range,
+        deuterium_flow,
+        levels=15,
+        cmap="seismic",
+        norm=plt.matplotlib.colors.CenteredNorm(vcenter=0),
+    )
+    axis.contour(
+        fuelling_range,
+        recycling_range,
+        deuterium_flow,
+        levels=[0],
+        colors="black",
+        linewidths=2,
+    )
+
+    # Plot star for mfile values
+    recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan)
+    fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan)
+    axis.plot(
+        fuelling_mfile,
+        recycling_mfile,
+        marker="*",
+        markersize=15,
+        color="yellow",
+        markeredgecolor="black",
+        markeredgewidth=1.5,
+    )
+
+    axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)")
+    axis.set_title("Plasma Deuterium Flow Rate (particles/s)", pad=20)
+    axis.minorticks_on()
+    axis.grid(True, which="major", linestyle="-", alpha=0.7)
+    axis.grid(True, which="minor", linestyle=":", alpha=0.4)
+    cbar = plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate")
+    cbar.ax.axhline(y=0, color="black", linewidth=2)
+
+
+def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
+    """Plot contour of alpha particle flow rate vs recycling and fuelling rate."""
+    fusion_dt_range = np.linspace(1e19, 5e21, 20)
+    f_t_alpha_energy_confinement_range = np.linspace(2.0, 10.0, 20)
+    alpha_flow = np.zeros((
+        len(fusion_dt_range),
+        len(f_t_alpha_energy_confinement_range),
+    ))
+
+    for i, fusion_dt in enumerate(fusion_dt_range):
+        for j, f_t_alpha_energy_confinement in enumerate(
+            f_t_alpha_energy_confinement_range
+        ):
+            alpha_flow[i, j] = PlasmaFuelling.calculate_plasma_alphas_flow_rate(
+                fusrat_dt_total=fusion_dt,
+                fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
+                t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
+                nd_plasma_alphas_vol_avg=mfile.get(
+                    "nd_plasma_alphas_vol_avg", scan=scan
+                ),
+                vol_plasma=mfile.get("vol_plasma", scan=scan),
+                f_t_alpha_energy_confinement=f_t_alpha_energy_confinement,
+            )
+
+    contour = axis.contourf(
+        f_t_alpha_energy_confinement_range,
+        fusion_dt_range,
+        alpha_flow,
+        levels=15,
+        cmap="seismic",
+        norm=plt.matplotlib.colors.CenteredNorm(vcenter=0),
+    )
+    axis.contour(
+        f_t_alpha_energy_confinement_range,
+        fusion_dt_range,
+        alpha_flow,
+        levels=[0],
+        colors="black",
+        linewidths=2,
+    )
+
+    # Plot star for mfile values
+    fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan)
+    f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan)
+    axis.plot(
+        f_t_alpha_mfile,
+        fusion_dt_mfile,
+        marker="*",
+        markersize=15,
+        color="yellow",
+        markeredgecolor="black",
+        markeredgewidth=1.5,
+    )
+
+    axis.set_xlabel(
+        "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)"
+    )
+    axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]")
+    axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)", pad=20)
+    axis.minorticks_on()
+    axis.grid(True, which="major", linestyle="-", alpha=0.7)
+    axis.grid(True, which="minor", linestyle=":", alpha=0.4)
+    cbar = plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate")
+    cbar.ax.axhline(y=0, color="black", linewidth=2)
+
+
+def plot_helium3_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
+    """Plot contour of helium-3 flow rate vs recycling and fuelling rate."""
+    recycling_range = np.linspace(0.01, 0.99, 20)
+    fuelling_range = np.linspace(0.01, 1.0, 20)
+    helium3_flow = np.zeros((len(recycling_range), len(fuelling_range)))
+
+    for i, recycling in enumerate(recycling_range):
+        for j, fuelling in enumerate(fuelling_range):
+            helium3_flow[i, j] = PlasmaFuelling.calculate_plasma_helium3_flow_rate(
+                f_molflow_plasma_fuelling_helium3=mfile.get(
+                    "f_molflow_plasma_fuelling_helium3", scan=scan
+                ),
+                eta_plasma_fuelling=fuelling,
+                molflow_plasma_fuelling_vv_injected=mfile.get(
+                    "molflow_plasma_fuelling_vv_injected", scan=scan
+                ),
+                fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
+                t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
+                f_plasma_particles_lcfs_recycled=recycling,
+                nd_plasma_fuel_ions_vol_avg=mfile.get(
+                    "nd_plasma_fuel_ions_vol_avg", scan=scan
+                ),
+                vol_plasma=mfile.get("vol_plasma", scan=scan),
+                f_plasma_fuel_helium3=mfile.get("f_plasma_fuel_helium3", scan=scan),
+            )
+
+    contour = axis.contourf(
+        fuelling_range,
+        recycling_range,
+        helium3_flow,
+        levels=15,
+        cmap="seismic",
+        norm=plt.matplotlib.colors.CenteredNorm(vcenter=0),
+    )
+    axis.contour(
+        fuelling_range,
+        recycling_range,
+        helium3_flow,
+        levels=[0],
+        colors="black",
+        linewidths=2,
+    )
+
+    # Plot star for mfile values
+    recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan)
+    fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan)
+    axis.plot(
+        fuelling_mfile,
+        recycling_mfile,
+        marker="*",
+        markersize=15,
+        color="yellow",
+        markeredgecolor="black",
+    )
+    axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)")
+    cbar = plt.colorbar(contour, ax=axis, label="Helium-3 Flow Rate")
+    cbar.ax.axhline(y=0, color="black", linewidth=2)
+    axis.set_title("Plasma Helium-3 Flow Rate (particles/s)", pad=20)
+    axis.minorticks_on()
+    axis.grid(True, which="major", linestyle="-", alpha=0.7)
+    axis.grid(True, which="minor", linestyle=":", alpha=0.4)
+
+
+def plot_fuelling_info(fig: plt.Figure, mfile: MFile, scan: int):
+    """Plot fuelling information."""
+    msg = (
+        f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n"
+        f"Total fuelling rate:"
+        f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n"
+        f"Total fuelling rate: "
+        f"{mfile.get('molflow_plasma_fuelling_vv_injected_moles', scan=scan):.4e} moles/s\n"
+        f"Total fuelling loss: "
+        f"{mfile.get('molflow_plasma_fuelling_loss', scan=scan):.4e} particles/s\n"
+        f"Total fuelling loss: "
+        f"{mfile.get('molflow_plasma_fuelling_loss_moles', scan=scan):.4e} moles/s\n"
+        f"Fuelling Rate Efficiency ($\\eta_{{\\text{{fuelling}}}}$): "
+        f"{mfile.get('eta_plasma_fuelling', scan=scan):.4f}\n"
+        f"Recycling Fraction ($R$): "
+        f"{mfile.get('f_plasma_particles_lcfs_recycled', scan=scan):.4f}\n\n"
+        f"Fraction of Tritium Fuelling: "
+        f"{mfile.get('f_molflow_plasma_fuelling_tritium', scan=scan):.4f}\n"
+        f"Fraction of Deuterium Fuelling: "
+        f"{mfile.get('f_molflow_plasma_fuelling_deuterium', scan=scan):.4f}\n"
+        f"Fraction of 3-Helium Fuelling: "
+        f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}\n\n"
+        f"Total Fuel Burnup Fraction: "
+        f"{mfile.get('f_plasma_fuel_burnup', scan=scan):.4f}\n"
+        f"Tritium Burnup Fraction: "
+        f"{mfile.get('f_plasma_tritium_burnup', scan=scan):.4f}\n"
+        f"Deuterium Burnup Fraction: "
+        f"{mfile.get('f_plasma_deuterium_burnup', scan=scan):.4f}"
+    )
+    fig.text(
+        0.75,
+        0.25,
+        msg,
+        ha="center",
+        va="center",
+        transform=fig.transFigure,
+        fontsize=9,
+        bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0},
+    )
+
+
 def main_plot(
     m_file: MFile,
     scan: int,
@@ -16647,6 +16961,21 @@ def _add_page(name: str | None = None):
     if m_file.get("i_plasma_shape", scan=scan) == PlasmaShapeModelType.SAUTER:
         plot_fusion_rate_contours(pages["rx_1_2"], pages["rx_3_4"], m_file, scan)
 
+    plot_tritium_flow_contour(
+        axis=_add_page("fuelling").add_subplot(231), mfile=m_file, scan=scan
+    )
+    plot_deuterium_flow_contour(
+        axis=pages["fuelling"].add_subplot(232), mfile=m_file, scan=scan
+    )
+    plot_helium3_flow_contour(
+        axis=pages["fuelling"].add_subplot(233), mfile=m_file, scan=scan
+    )
+    plot_alpha_flow_contour(
+        axis=pages["fuelling"].add_subplot(223), mfile=m_file, scan=scan
+    )
+
+    plot_fuelling_info(pages["fuelling"], m_file, scan)
+
     if i_shape != PlasmaShapeModelType.SAUTER:
         msg = (
             "Fusion-rate contour plots require a closed (Sauter) plasma boundary "

From d8f687265acd30558f505e450815fad07a700bdd Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 13:19:54 +0100
Subject: [PATCH 16/39] Add new constraints to large_tokamak input and fix some
 runtime bugs

---
 process/core/io/plot/summary.py               |  4 +-
 process/core/solver/constraints.py            |  2 +-
 process/main.py                               |  2 +-
 process/models/physics/fuelling.py            |  3 +-
 .../input_files/large_tokamak_nof.IN.DAT      | 38 ++++++++++++++++++-
 5 files changed, 42 insertions(+), 7 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index f903cd97b0..41c4bd471f 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -12228,7 +12228,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
 
     axis.text(
         0.285,
-        0.56,
+        0.58,
         "$\\text{D - 3He}$",
         fontsize=20,
         verticalalignment="top",
@@ -16661,7 +16661,7 @@ def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
 
     # Plot star for mfile values
     fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan)
-    f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan)
+    f_t_alpha_mfile = mfile.get("f_t_alpha_energy_confinement", scan=scan)
     axis.plot(
         f_t_alpha_mfile,
         fusion_dt_mfile,
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 5b68c9194a..1dc21982f8 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2079,7 +2079,7 @@ def constraint_equation_96(constraint_registration, data):
     # Alpha particle balance
     numerator = data.physics.fusrat_dt_total + data.physics.fusrat_plasma_dhe3
     denominator = (data.physics.nd_plasma_alphas_vol_avg * data.physics.vol_plasma) / (
-        data.physics.t_energy_confinement * data.physics.f_alpha_energy_confinement
+        data.physics.t_energy_confinement * data.physics.f_t_alpha_energy_confinement
     )
 
     return eq(numerator, denominator, constraint_registration)
diff --git a/process/main.py b/process/main.py
index 4e282de5ff..356a23958a 100644
--- a/process/main.py
+++ b/process/main.py
@@ -827,7 +827,7 @@ def models(self) -> tuple[Model, ...]:
             self.plasma_fields,
             self.sauter_bootstrap_current,
             self.plasma_transition,
-            self.self.plasma_fuelling,
+            self.plasma_fuelling,
             self.physics_detailed,
             self.electron_cyclotron,
             self.lower_hybrid,
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index b9e70799b8..fd8bc1c2e7 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -2,7 +2,6 @@
 from process.core import process_output as po
 from process.core.model import Model
 from process.data_structure import (
-    numerics,
     reinke_variables,
 )
 
@@ -459,7 +458,7 @@ def output_fuelling_info(self):
             "OP ",
         )
 
-        if 78 in numerics.icc:
+        if 78 in self.data.numerics.icc:
             po.osubhd(self.outfile, "Reinke Criterion :")
             po.ovarin(
                 self.outfile,
diff --git a/tests/regression/input_files/large_tokamak_nof.IN.DAT b/tests/regression/input_files/large_tokamak_nof.IN.DAT
index 828a3a6e34..7cfd5e5ced 100644
--- a/tests/regression/input_files/large_tokamak_nof.IN.DAT
+++ b/tests/regression/input_files/large_tokamak_nof.IN.DAT
@@ -19,7 +19,7 @@ epsvmc = 1e-7
 
 * Number of equality constraints
 * (the first n_equality_constraints icc=... statements are equality constraints)
-n_equality_constraints = 3
+n_equality_constraints = 7
 
 * Maximum number of solver iterations
 maxcal = 100
@@ -41,6 +41,18 @@ icc = 2
 *--------------------------*
 icc = 11
 
+* Tritium particle balance
+icc = 93
+
+* Deuterium particle balance
+icc = 94
+
+* Alpha particle balance
+icc = 96
+
+* Fuelling composition consistency
+icc = 97
+
 * Constraint Equations - Limit Equations *
 ******************************************
 
@@ -178,6 +190,30 @@ t_burn_min = 7200.0
 * Iteration Variables *
 ***********************
 
+* Particle recycling fraction
+ixc = 177
+f_plasma_particles_lcfs_recycled = 0.9
+
+* Plasma fuelling efficiecy
+ixc = 178
+eta_plasma_fuelling = 0.7
+
+* Injected VV fuelling rate
+ixc = 179
+molflow_plasma_fuelling_vv_injected = 5e21
+boundl(179) = 1e20
+
+* Deuterium fuelling fraction
+ixc = 180
+f_molflow_plasma_fuelling_deuterium = 0.5
+boundl(180) = 0.4
+
+* Tritium fuelling fraction
+ixc = 181
+f_molflow_plasma_fuelling_tritium = 0.5
+boundl(181) = 0.4
+
+
 * b_plasma_toroidal_on_axis [T]
 ixc = 2
 b_plasma_toroidal_on_axis = 5.7

From 41ea5cb031ecc966416084847812e3a10c7febd2 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 1 Jul 2026 08:49:43 +0100
Subject: [PATCH 17/39] Add more descriptive notes to each function about what
 the output represents

---
 process/core/io/plot/summary.py    |  4 +--
 process/models/physics/fuelling.py | 54 ++++++++++++++++++++----------
 2 files changed, 39 insertions(+), 19 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 41c4bd471f..260f095f12 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -12241,8 +12241,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
         f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.4f} MW\n"
         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"
+        f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m$^3$/s\n"
+        f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m$^3$/s\n\n"
         f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m$^3$\n"
         f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m$^3$\n\n"
         f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m$^3$\n"
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index fd8bc1c2e7..4bfac0c5e2 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -7,7 +7,7 @@
 
 
 class PlasmaFuelling(Model):
-    """Class to hold plasma fuelling calculations for plasma processing."""
+    """Class to hold plasma fuelling calculations and output."""
 
     def __init__(self):
 
@@ -157,7 +157,7 @@ def calculate_plasma_tritium_flow_rate(
         vol_plasma: float,
         f_plasma_fuel_tritium: float,
     ) -> float:
-        """Calculate the tritium flow rate in the plasma exhaust.
+        """Calculate the tritium flow rate into the plasma.
 
         Parameters
         ----------
@@ -176,16 +176,21 @@ def calculate_plasma_tritium_flow_rate(
         f_plasma_particles_lcfs_recycled : float
             Fraction of plasma particles recycled at the LCFS.
         nd_plasma_fuel_ions_vol_avg : float
-            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
         vol_plasma : float
-            Plasma volume (m^3).
+            Plasma volume (m³).
         f_plasma_fuel_tritium : float
             Fraction of tritium in the plasma fuel.
 
         Returns
         -------
         float
-            Tritium flow rate in the plasma exhaust (particles/s).
+            Tritium flow rate into the plasma (particles/s).
+
+        Notes
+        -----
+        - A positive value indicates a net flow of tritium into the plasma,
+        while a negative value indicates a net loss of tritium from the plasma.
 
         """
         return (
@@ -216,7 +221,7 @@ def calculate_plasma_deuterium_flow_rate(
         vol_plasma: float,
         f_plasma_fuel_deuterium: float,
     ) -> float:
-        """Calculate the deuterium flow rate in the plasma exhaust.
+        """Calculate the deuterium flow rate into the plasma.
 
         Parameters
         ----------
@@ -237,16 +242,21 @@ def calculate_plasma_deuterium_flow_rate(
         f_plasma_particles_lcfs_recycled : float
             Fraction of plasma particles recycled at the LCFS.
         nd_plasma_fuel_ions_vol_avg : float
-            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
         vol_plasma : float
-            Plasma volume (m^3).
+            Plasma volume (m³).
         f_plasma_fuel_deuterium : float
             Fraction of deuterium in the plasma fuel.
 
         Returns
         -------
         float
-            Deuterium flow rate in the plasma exhaust (particles/s).
+            Deuterium flow rate into the plasma (particles/s).
+
+        Notes
+        -----
+        - A positive value indicates a net flow of deuterium into the plasma,
+        while a negative value indicates a net loss of deuterium from the plasma.
 
 
         """
@@ -277,7 +287,7 @@ def calculate_plasma_helium3_flow_rate(
         vol_plasma: float,
         f_plasma_fuel_helium3: float,
     ) -> float:
-        """Calculate the helium-3 flow rate in the plasma exhaust.
+        """Calculate the helium-3 flow rate into the plasma.
 
         Parameters
         ----------
@@ -294,16 +304,21 @@ def calculate_plasma_helium3_flow_rate(
         f_plasma_particles_lcfs_recycled : float
             Fraction of plasma particles recycled at the LCFS.
         nd_plasma_fuel_ions_vol_avg : float
-            Volume-averaged density of fuel ions in the plasma (particles/m^3).
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
         vol_plasma : float
-            Plasma volume (m^3).
+            Plasma volume (m³).
         f_plasma_fuel_helium3 : float
             Fraction of helium-3 in the plasma fuel.
 
         Returns
         -------
         float
-            Helium-3 flow rate in the plasma exhaust (particles/s).
+            Helium-3 flow rate into the plasma (particles/s).
+
+        Notes
+        -----
+        - A positive value indicates a net flow of helium-3 into the plasma,
+        while a negative value indicates a net loss of helium-3 from the plasma.
 
         """
         return (
@@ -328,7 +343,7 @@ def calculate_plasma_alphas_flow_rate(
         nd_plasma_alphas_vol_avg: float,
         vol_plasma: float,
     ) -> float:
-        """Calculate the alpha particle flow rate in the plasma exhaust.
+        """Calculate the net alpha particle flow rate into the plasma.
 
         Parameters
         ----------
@@ -341,14 +356,19 @@ def calculate_plasma_alphas_flow_rate(
         f_t_alpha_energy_confinement : float
             Ratio of alpha particle confinement time to energy confinement time (dimensionless).
         nd_plasma_alphas_vol_avg : float
-            Volume-averaged density of alpha particles in the plasma (particles/m^3).
+            Volume-averaged density of alpha particles in the plasma (particles/m³).
         vol_plasma : float
-            Plasma volume (m^3).
+            Plasma volume (m³).
 
         Returns
         -------
         float
-            Alpha particle flow rate in the plasma exhaust (particles/s).
+            Alpha particle flow rate into the plasma (particles/s).
+
+        Notes
+        -----
+        - A positive value indicates a net flow of alpha particles into the plasma,
+        while a negative value indicates a net loss of alpha particles from the plasma.
 
         """
         # Alpha particle balance

From b1cfccd29ba82efa33842f6db62c64c692b7bf30 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 1 Jul 2026 16:11:55 +0100
Subject: [PATCH 18/39] Refactor plasma fuelling calculations to use dedicated
 methods for tritium source and loss rates. Implement these new methods in the
 constraints

---
 process/core/solver/constraints.py | 30 +++++-----
 process/models/physics/fuelling.py | 95 ++++++++++++++++++++++++++----
 2 files changed, 97 insertions(+), 28 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 1dc21982f8..9df7e7a66d 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -14,6 +14,7 @@
 from process.data_structure.build_variables import TFCSRadialConfiguration
 from process.models.physics.density_limit import DensityLimitModel
 from process.models.physics.exhaust import PlasmaExhaust
+from process.models.physics.fuelling import PlasmaFuelling
 from process.models.physics.physics import (
     BetaComponentLimits,
     PlasmaBeta,
@@ -1970,25 +1971,22 @@ def constraint_equation_93(constraint_registration, data):
     The equation ensures that the rate of tritium production equals the rate of tritium
     loss, maintaining a steady-state condition for tritium in the plasma.
     """
-    numerator = (
-        (
-            data.physics.eta_plasma_fuelling
-            * data.physics.molflow_plasma_fuelling_vv_injected
-            * data.physics.f_molflow_plasma_fuelling_tritium
-        )
-        + data.physics.fusrat_plasma_dd_triton
-        - data.physics.fusrat_dt_total
+    source = PlasmaFuelling.calculate_plasma_tritium_source_rate(
+        f_molflow_plasma_fuelling_tritium=data.physics.f_molflow_plasma_fuelling_tritium,
+        eta_plasma_fuelling=data.physics.eta_plasma_fuelling,
+        molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected,
+        fusrat_plasma_dd_triton=data.physics.fusrat_plasma_dd_triton,
     )
-    denominator = (
-        data.physics.nd_plasma_fuel_ions_vol_avg
-        * data.physics.vol_plasma
-        * data.physics.f_plasma_fuel_tritium
-    ) / (
-        data.physics.t_energy_confinement
-        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    sink = PlasmaFuelling.calculate_plasma_tritium_loss_rate(
+        fusrat_dt_total=data.physics.fusrat_dt_total,
+        t_energy_confinement=data.physics.t_energy_confinement,
+        f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled,
+        nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg,
+        vol_plasma=data.physics.vol_plasma,
+        f_plasma_fuel_tritium=data.physics.f_plasma_fuel_tritium,
     )
 
-    return eq(numerator, denominator, constraint_registration)
+    return eq(source, -sink, constraint_registration)
 
 
 @ConstraintManager.register_constraint(94, "particles/s", "=")
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 4bfac0c5e2..64317ae0e1 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -144,8 +144,8 @@ def calculate_deuterium_burnup_fraction(
             molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium
         )
 
-    @staticmethod
     def calculate_plasma_tritium_flow_rate(
+        self,
         f_molflow_plasma_fuelling_tritium: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
@@ -192,19 +192,90 @@ def calculate_plasma_tritium_flow_rate(
         - A positive value indicates a net flow of tritium into the plasma,
         while a negative value indicates a net loss of tritium from the plasma.
 
+        """
+        return self.calculate_plasma_tritium_source_rate(
+            f_molflow_plasma_fuelling_tritium=f_molflow_plasma_fuelling_tritium,
+            eta_plasma_fuelling=eta_plasma_fuelling,
+            molflow_plasma_fuelling_vv_injected=molflow_plasma_fuelling_vv_injected,
+            fusrat_plasma_dd_triton=fusrat_plasma_dd_triton,
+        ) + self.calculate_plasma_tritium_loss_rate(
+            fusrat_dt_total=fusrat_dt_total,
+            t_energy_confinement=t_energy_confinement,
+            f_plasma_particles_lcfs_recycled=f_plasma_particles_lcfs_recycled,
+            nd_plasma_fuel_ions_vol_avg=nd_plasma_fuel_ions_vol_avg,
+            vol_plasma=vol_plasma,
+            f_plasma_fuel_tritium=f_plasma_fuel_tritium,
+        )
+
+    @staticmethod
+    def calculate_plasma_tritium_source_rate(
+        f_molflow_plasma_fuelling_tritium: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        fusrat_plasma_dd_triton: float,
+    ) -> float:
+        """Calculate the tritium source rate in the plasma.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_tritium : float
+            Fraction of tritium in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+        fusrat_plasma_dd_triton : float
+            Tritium production rate from D-D fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+
+        Returns
+        -------
+        float
+            Tritium source rate in the plasma (particles/s).
+
         """
         return (
-            (
-                f_molflow_plasma_fuelling_tritium
-                * eta_plasma_fuelling
-                * molflow_plasma_fuelling_vv_injected
-            )
-            - fusrat_dt_total
-            + fusrat_plasma_dd_triton
-            - (
-                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium)
-                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
-            )
+            f_molflow_plasma_fuelling_tritium
+            * eta_plasma_fuelling
+            * molflow_plasma_fuelling_vv_injected
+        ) + fusrat_plasma_dd_triton
+
+    @staticmethod
+    def calculate_plasma_tritium_loss_rate(
+        fusrat_dt_total: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_tritium: float,
+    ) -> float:
+        """Calculate the tritium loss rate from the plasma.
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
+        vol_plasma : float
+            Plasma volume (m³).
+        f_plasma_fuel_tritium : float
+            Fraction of tritium in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Tritium loss rate from the plasma (particles/s).
+
+        """
+        return -fusrat_dt_total - (
+            (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium)
+            / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
         )
 
     @staticmethod

From b92a538d5a593bed40987d41f756db04848978e8 Mon Sep 17 00:00:00 2001
From: Christopher Ashe <91618944+chris-ashe@users.noreply.github.com>
Date: Wed, 22 Jul 2026 11:41:33 +0100
Subject: [PATCH 19/39] Update
 documentation/source/physics-models/plasma_fuelling.md

Co-authored-by: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com>
---
 documentation/source/physics-models/plasma_fuelling.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index b9a05389e3..0521aa6b3e 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -26,7 +26,7 @@ $$
 \frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha}}{dt} = 0
 $$
 
-Here $\eta_{\text{fuelling}}$ is the fuelling efficiecny which represents the method of injecting fuel into the plasma. Gas puffing on the low field side is probably around 0.01-0.1, supersonic gas is 0.1 and 0.2 and using pellets can get you close to unity with 0.5-0.9. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel, so $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ together presents the fraction of injected fuel that actually makes it into the plasma core to fuse. 
+Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions.
 
 
 The fuelling fractional compositions is given by $f$

From e9fc55e37bd037ef096c5acd11e7d3c506c7bd88 Mon Sep 17 00:00:00 2001
From: Christopher Ashe <91618944+chris-ashe@users.noreply.github.com>
Date: Wed, 22 Jul 2026 15:27:57 +0100
Subject: [PATCH 20/39] Update
 documentation/source/physics-models/plasma_fuelling.md

Co-authored-by: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com>
---
 documentation/source/physics-models/plasma_fuelling.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index 0521aa6b3e..e3f251b4fe 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -109,7 +109,7 @@ $$
 f_b = \frac{\Gamma_{\text{fusion}}}{\Gamma_{\text{fuel}}}
 $$
 
-where $\Gamma_{\text{fusion}}$ is the fusion rate in reactions per second $[\text{s}^{-1}]$, and $\Gamma_{\text{fuel}}$ is our fuelling rate into the vessel like is shown above also in $[\text{s}^{-1}]$.
+where $\Gamma_{\text{fusion}}$ is the fusion reaction rate per second $[\text{s}^{-1}]$, and $\Gamma_{\text{fuel}}$ is the fuel injection rate into the vacuum vessel in particles per second $[\text{s}^{-1}]$.
 
 -----------------
 

From 484fea1b1114b40ce1e5ca5ce8a823f1d0b63284 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 22 Jul 2026 16:00:55 +0100
Subject: [PATCH 21/39] Quick fixes to correct for new thermal naming of alphas

---
 process/core/io/plot/summary.py    | 7 ++++---
 process/core/solver/constraints.py | 6 +++---
 process/models/physics/fuelling.py | 8 +++-----
 process/models/physics/physics.py  | 2 +-
 4 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 260f095f12..a9826e9013 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16491,7 +16491,8 @@ def plot_tritium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
 
     for i, recycling in enumerate(recycling_range):
         for j, fuelling in enumerate(fuelling_range):
-            tritium_flow[i, j] = PlasmaFuelling.calculate_plasma_tritium_flow_rate(
+            plasma_fuelling = PlasmaFuelling()
+            tritium_flow[i, j] = plasma_fuelling.calculate_plasma_tritium_flow_rate(
                 f_molflow_plasma_fuelling_tritium=mfile.get(
                     "f_molflow_plasma_fuelling_tritium", scan=scan
                 ),
@@ -16635,8 +16636,8 @@ def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
                 fusrat_dt_total=fusion_dt,
                 fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
                 t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
-                nd_plasma_alphas_vol_avg=mfile.get(
-                    "nd_plasma_alphas_vol_avg", scan=scan
+                nd_plasma_alphas_thermal_vol_avg=mfile.get(
+                    "nd_plasma_alphas_thermal_vol_avg", scan=scan
                 ),
                 vol_plasma=mfile.get("vol_plasma", scan=scan),
                 f_t_alpha_energy_confinement=f_t_alpha_energy_confinement,
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 9df7e7a66d..4bcdcc5e77 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2076,9 +2076,9 @@ def constraint_equation_96(constraint_registration, data):
     """
     # Alpha particle balance
     numerator = data.physics.fusrat_dt_total + data.physics.fusrat_plasma_dhe3
-    denominator = (data.physics.nd_plasma_alphas_vol_avg * data.physics.vol_plasma) / (
-        data.physics.t_energy_confinement * data.physics.f_t_alpha_energy_confinement
-    )
+    denominator = (
+        data.physics.nd_plasma_alphas_thermal_vol_avg * data.physics.vol_plasma
+    ) / (data.physics.t_energy_confinement * data.physics.f_t_alpha_energy_confinement)
 
     return eq(numerator, denominator, constraint_registration)
 
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 64317ae0e1..b0998baf98 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -226,8 +226,6 @@ def calculate_plasma_tritium_source_rate(
             Total fuelling rate (particles/s).
         fusrat_plasma_dd_triton : float
             Tritium production rate from D-D fusion (particles/s).
-        t_energy_confinement : float
-            Energy confinement time (s).
 
         Returns
         -------
@@ -411,7 +409,7 @@ def calculate_plasma_alphas_flow_rate(
         fusrat_plasma_dhe3: float,
         t_energy_confinement: float,
         f_t_alpha_energy_confinement: float,
-        nd_plasma_alphas_vol_avg: float,
+        nd_plasma_alphas_thermal_vol_avg: float,
         vol_plasma: float,
     ) -> float:
         """Calculate the net alpha particle flow rate into the plasma.
@@ -426,7 +424,7 @@ def calculate_plasma_alphas_flow_rate(
             Energy confinement time (s).
         f_t_alpha_energy_confinement : float
             Ratio of alpha particle confinement time to energy confinement time (dimensionless).
-        nd_plasma_alphas_vol_avg : float
+        nd_plasma_alphas_thermal_vol_avg : float
             Volume-averaged density of alpha particles in the plasma (particles/m³).
         vol_plasma : float
             Plasma volume (m³).
@@ -447,7 +445,7 @@ def calculate_plasma_alphas_flow_rate(
         return (
             fusrat_dt_total
             + fusrat_plasma_dhe3
-            - (nd_plasma_alphas_vol_avg * vol_plasma)
+            - (nd_plasma_alphas_thermal_vol_avg * vol_plasma)
             / (t_energy_confinement * f_t_alpha_energy_confinement)
         )
 
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 9b2fb36588..79ef1f58cc 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -979,7 +979,7 @@ def run(self):
             self.data.physics.f_t_alpha_energy_confinement,
         ) = self.phyaux(
             fusden_alpha_total=self.data.physics.fusden_alpha_total,
-            nd_plasma_alphas_vol_avg=self.data.physics.nd_plasma_alphas_thermal_vol_avg,
+            nd_plasma_alphas_thermal_vol_avg=self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             t_energy_confinement=self.data.physics.t_energy_confinement,
         )
 

From ca9da06a956bb121f90b99a13ca8efd722c44c8a Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 08:54:11 +0100
Subject: [PATCH 22/39] Refactor plasma fuelling calculations to use dedicated
 methods for deuterium source and loss rates

---
 process/core/io/plot/summary.py    |   2 +-
 process/core/solver/constraints.py |  31 ++++-----
 process/models/physics/fuelling.py | 103 +++++++++++++++++++++++++----
 3 files changed, 105 insertions(+), 31 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index a9826e9013..3f042f205a 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16560,7 +16560,7 @@ def plot_deuterium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
 
     for i, recycling in enumerate(recycling_range):
         for j, fuelling in enumerate(fuelling_range):
-            deuterium_flow[i, j] = PlasmaFuelling.calculate_plasma_deuterium_flow_rate(
+            deuterium_flow[i, j] = PlasmaFuelling().calculate_plasma_deuterium_flow_rate(
                 f_molflow_plasma_fuelling_deuterium=mfile.get(
                     "f_molflow_plasma_fuelling_deuterium", scan=scan
                 ),
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 4bcdcc5e77..2116190b5e 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2004,26 +2004,23 @@ def constraint_equation_94(constraint_registration, data):
     The equation ensures that the rate of deuterium production equals the rate of deuterium
     loss, maintaining a steady-state condition for deuterium in the plasma.
     """
-    numerator = (
-        (
-            data.physics.eta_plasma_fuelling
-            * data.physics.molflow_plasma_fuelling_vv_injected
-            * data.physics.f_molflow_plasma_fuelling_deuterium
-        )
-        - data.physics.fusrat_dt_total
-        - data.physics.fusrat_plasma_dhe3
-        - 2.0 * data.physics.fusrat_plasma_dd_total
+    source = PlasmaFuelling.calculate_plasma_deuterium_source_rate(
+        f_molflow_plasma_fuelling_deuterium=data.physics.f_molflow_plasma_fuelling_deuterium,
+        eta_plasma_fuelling=data.physics.eta_plasma_fuelling,
+        molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected,
     )
-    denominator = (
-        data.physics.nd_plasma_fuel_ions_vol_avg
-        * data.physics.vol_plasma
-        * data.physics.f_plasma_fuel_deuterium
-    ) / (
-        data.physics.t_energy_confinement
-        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    sink = PlasmaFuelling.calculate_plasma_deuterium_loss_rate(
+        fusrat_dt_total=data.physics.fusrat_dt_total,
+        fusrat_plasma_dd_total=data.physics.fusrat_plasma_dd_total,
+        fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3,
+        t_energy_confinement=data.physics.t_energy_confinement,
+        f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled,
+        nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg,
+        vol_plasma=data.physics.vol_plasma,
+        f_plasma_fuel_deuterium=data.physics.f_plasma_fuel_deuterium,
     )
 
-    return eq(numerator, denominator, constraint_registration)
+    return eq(source, -sink, constraint_registration)
 
 
 @ConstraintManager.register_constraint(95, "particles/s", "=")
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index b0998baf98..18b7d978d0 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -277,7 +277,84 @@ def calculate_plasma_tritium_loss_rate(
         )
 
     @staticmethod
+    def calculate_plasma_deuterium_source_rate(
+        f_molflow_plasma_fuelling_deuterium: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+    ) -> float:
+        """Calculate the deuterium source rate in the plasma.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_deuterium : float
+            Fraction of deuterium in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+
+        Returns
+        -------
+        float
+            Deuterium source rate in the plasma (particles/s).
+
+        """
+        return (
+            f_molflow_plasma_fuelling_deuterium
+            * eta_plasma_fuelling
+            * molflow_plasma_fuelling_vv_injected
+        )
+
+    @staticmethod
+    def calculate_plasma_deuterium_loss_rate(
+        fusrat_dt_total: float,
+        fusrat_plasma_dd_total: float,
+        fusrat_plasma_dhe3: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_deuterium: float,
+    ) -> float:
+        """Calculate the deuterium loss rate from the plasma.
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        fusrat_plasma_dd_total : float
+            Total deuterium consumption rate from DD fusion (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
+        vol_plasma : float
+            Plasma volume (m³).
+        f_plasma_fuel_deuterium : float
+            Fraction of deuterium in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Deuterium loss rate from the plasma (particles/s).
+
+        """
+        return (
+            -fusrat_dt_total
+            - 2 * fusrat_plasma_dd_total
+            - fusrat_plasma_dhe3
+            - (
+                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium)
+                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
+            )
+        )
+
     def calculate_plasma_deuterium_flow_rate(
+        self,
         f_molflow_plasma_fuelling_deuterium: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
@@ -329,19 +406,19 @@ def calculate_plasma_deuterium_flow_rate(
 
 
         """
-        return (
-            (
-                f_molflow_plasma_fuelling_deuterium
-                * eta_plasma_fuelling
-                * molflow_plasma_fuelling_vv_injected
-            )
-            - fusrat_dt_total
-            - 2 * fusrat_plasma_dd_total
-            - fusrat_plasma_dhe3
-            - (
-                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium)
-                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
-            )
+        return self.calculate_plasma_deuterium_source_rate(
+            f_molflow_plasma_fuelling_deuterium=f_molflow_plasma_fuelling_deuterium,
+            eta_plasma_fuelling=eta_plasma_fuelling,
+            molflow_plasma_fuelling_vv_injected=molflow_plasma_fuelling_vv_injected,
+        ) + self.calculate_plasma_deuterium_loss_rate(
+            fusrat_dt_total=fusrat_dt_total,
+            fusrat_plasma_dd_total=fusrat_plasma_dd_total,
+            fusrat_plasma_dhe3=fusrat_plasma_dhe3,
+            t_energy_confinement=t_energy_confinement,
+            f_plasma_particles_lcfs_recycled=f_plasma_particles_lcfs_recycled,
+            nd_plasma_fuel_ions_vol_avg=nd_plasma_fuel_ions_vol_avg,
+            vol_plasma=vol_plasma,
+            f_plasma_fuel_deuterium=f_plasma_fuel_deuterium,
         )
 
     @staticmethod

From af78c83be5410d42c7f779ed3f1a22d686fa38d4 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 09:33:18 +0100
Subject: [PATCH 23/39] Implement helium-3 source and loss rate calculations in
 PlasmaFuelling class; update constraints and flow rate calculations
 accordingly.

---
 .../source/physics-models/plasma_fuelling.md  |  4 +-
 process/core/io/plot/summary.py               |  3 +-
 process/core/solver/constraints.py            | 27 +++---
 process/models/physics/fuelling.py            | 96 ++++++++++++++++---
 4 files changed, 103 insertions(+), 27 deletions(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index e3f251b4fe..f7c88a632c 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -13,7 +13,7 @@ $$
 $$
 
 $$
-\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} -\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
 $$
 
 $$
@@ -88,7 +88,7 @@ $$
 ### Helium-3 Flow Rate | `calculate_plasma_helium3_flow_rate()`
 
 $$
-\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}}-\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
 $$
 
 ---------------
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 3f042f205a..7112e50496 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16693,7 +16693,7 @@ def plot_helium3_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
 
     for i, recycling in enumerate(recycling_range):
         for j, fuelling in enumerate(fuelling_range):
-            helium3_flow[i, j] = PlasmaFuelling.calculate_plasma_helium3_flow_rate(
+            helium3_flow[i, j] = PlasmaFuelling().calculate_plasma_helium3_flow_rate(
                 f_molflow_plasma_fuelling_helium3=mfile.get(
                     "f_molflow_plasma_fuelling_helium3", scan=scan
                 ),
@@ -16702,6 +16702,7 @@ def plot_helium3_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
                     "molflow_plasma_fuelling_vv_injected", scan=scan
                 ),
                 fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
+                fusrat_plasma_dd_helion=mfile.get("fusrat_plasma_dd_helion", scan=scan),
                 t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
                 f_plasma_particles_lcfs_recycled=recycling,
                 nd_plasma_fuel_ions_vol_avg=mfile.get(
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 2116190b5e..d6be3e39e5 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2038,21 +2038,22 @@ def constraint_equation_95(constraint_registration, data):
     The equation ensures that the rate of helium-3 production equals the rate of helium-3
     loss, maintaining a steady-state condition for helium-3 in the plasma.
     """
-    numerator = (
-        data.physics.eta_plasma_fuelling
-        * data.physics.molflow_plasma_fuelling_vv_injected
-        * data.physics.f_molflow_plasma_fuelling_helium3
-    ) + data.physics.fusrat_plasma_dhe3
-    denominator = (
-        data.physics.nd_plasma_fuel_ions_vol_avg
-        * data.physics.vol_plasma
-        * data.physics.f_plasma_fuel_helium3
-    ) / (
-        data.physics.t_energy_confinement
-        / (1 - data.physics.f_plasma_particles_lcfs_recycled)
+    source = PlasmaFuelling.calculate_plasma_helium3_source_rate(
+        f_molflow_plasma_fuelling_helium3=data.physics.f_molflow_plasma_fuelling_helium3,
+        eta_plasma_fuelling=data.physics.eta_plasma_fuelling,
+        molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected,
+        fusrat_plasma_dd_helion=data.physics.fusrat_plasma_dd_helion,
+    )
+    sink = PlasmaFuelling.calculate_plasma_helium3_loss_rate(
+        fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3,
+        t_energy_confinement=data.physics.t_energy_confinement,
+        f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled,
+        nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg,
+        vol_plasma=data.physics.vol_plasma,
+        f_plasma_fuel_helium3=data.physics.f_plasma_fuel_helium3,
     )
 
-    return eq(numerator, denominator, constraint_registration)
+    return eq(source, -sink, constraint_registration)
 
 
 @ConstraintManager.register_constraint(96, "particles/s", "=")
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 18b7d978d0..111df2e745 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -422,11 +422,82 @@ def calculate_plasma_deuterium_flow_rate(
         )
 
     @staticmethod
+    def calculate_plasma_helium3_source_rate(
+        f_molflow_plasma_fuelling_helium3: float,
+        eta_plasma_fuelling: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        fusrat_plasma_dd_helion: float,
+    ) -> float:
+        """Calculate the helium-3 source rate in the plasma.
+
+        Parameters
+        ----------
+        f_molflow_plasma_fuelling_helium3 : float
+            Fraction of helium-3 in the plasma fuelling.
+        eta_plasma_fuelling : float
+            Fuelling rate efficiency.
+        molflow_plasma_fuelling_vv_injected : float
+            Total fuelling rate (particles/s).
+        fusrat_plasma_dd_helion : float
+            Helium-3 production rate from DD fusion (particles/s).
+
+        Returns
+        -------
+        float
+            Helium-3 source rate in the plasma (particles/s).
+
+        """
+        return (
+            f_molflow_plasma_fuelling_helium3
+            * eta_plasma_fuelling
+            * molflow_plasma_fuelling_vv_injected
+        ) + fusrat_plasma_dd_helion
+
+    @staticmethod
+    def calculate_plasma_helium3_loss_rate(
+        fusrat_plasma_dhe3: float,
+        t_energy_confinement: float,
+        f_plasma_particles_lcfs_recycled: float,
+        nd_plasma_fuel_ions_vol_avg: float,
+        vol_plasma: float,
+        f_plasma_fuel_helium3: float,
+    ) -> float:
+        """Calculate the helium-3 loss rate from the plasma.
+
+        Parameters
+        ----------
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_plasma_particles_lcfs_recycled : float
+            Fraction of plasma particles recycled at the LCFS.
+        nd_plasma_fuel_ions_vol_avg : float
+            Volume-averaged density of fuel ions in the plasma (particles/m³).
+        vol_plasma : float
+            Plasma volume (m³).
+        f_plasma_fuel_helium3 : float
+            Fraction of helium-3 in the plasma fuel.
+
+        Returns
+        -------
+        float
+            Helium-3 loss rate from the plasma (particles/s).
+
+
+        """
+        return -fusrat_plasma_dhe3 - (
+            (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_helium3)
+            / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
+        )
+
     def calculate_plasma_helium3_flow_rate(
+        self,
         f_molflow_plasma_fuelling_helium3: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
         fusrat_plasma_dhe3: float,
+        fusrat_plasma_dd_helion: float,
         t_energy_confinement: float,
         f_plasma_particles_lcfs_recycled: float,
         nd_plasma_fuel_ions_vol_avg: float,
@@ -445,6 +516,8 @@ def calculate_plasma_helium3_flow_rate(
             Total fuelling rate (particles/s).
         fusrat_plasma_dhe3 : float
             Deuterium consumption rate from D-He3 fusion (particles/s).
+        fusrat_plasma_dd_helion : float
+            Helium-3 production rate from DD fusion (particles/s).
         t_energy_confinement : float
             Energy confinement time (s).
         f_plasma_particles_lcfs_recycled : float
@@ -467,17 +540,18 @@ def calculate_plasma_helium3_flow_rate(
         while a negative value indicates a net loss of helium-3 from the plasma.
 
         """
-        return (
-            (
-                f_molflow_plasma_fuelling_helium3
-                * eta_plasma_fuelling
-                * molflow_plasma_fuelling_vv_injected
-            )
-            + fusrat_plasma_dhe3
-            - (
-                (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_helium3)
-                / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled))
-            )
+        return self.calculate_plasma_helium3_source_rate(
+            f_molflow_plasma_fuelling_helium3=f_molflow_plasma_fuelling_helium3,
+            eta_plasma_fuelling=eta_plasma_fuelling,
+            molflow_plasma_fuelling_vv_injected=molflow_plasma_fuelling_vv_injected,
+            fusrat_plasma_dd_helion=fusrat_plasma_dd_helion,
+        ) + self.calculate_plasma_helium3_loss_rate(
+            fusrat_plasma_dhe3=fusrat_plasma_dhe3,
+            t_energy_confinement=t_energy_confinement,
+            f_plasma_particles_lcfs_recycled=f_plasma_particles_lcfs_recycled,
+            nd_plasma_fuel_ions_vol_avg=nd_plasma_fuel_ions_vol_avg,
+            vol_plasma=vol_plasma,
+            f_plasma_fuel_helium3=f_plasma_fuel_helium3,
         )
 
     @staticmethod

From 8328cc9ec7c344ea6438d7daf6781663fbc29e11 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 09:45:15 +0100
Subject: [PATCH 24/39] Refactor alpha particle flow calculations to use
 thermal terminology; update related constraints and documentation
 accordingly.

---
 .../source/physics-models/plasma_fuelling.md  | 25 ++++++-------------
 process/core/io/plot/summary.py               |  2 +-
 process/core/solver/constraints.py            | 14 +++++------
 process/models/physics/fuelling.py            | 13 +++++-----
 4 files changed, 21 insertions(+), 33 deletions(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index f7c88a632c..3a49869334 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -2,7 +2,7 @@
 
 ## Particle balance
 
-The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$.
+The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the thermal helium $\alpha$ ash.
 
 $$
 \frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
@@ -17,13 +17,13 @@ $$
 $$
 
 $$
-\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*}
+\frac{dn_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
 $$
 
 In a steady state equilibrium all 4 of these equations should balance, therefore:
 
 $$
-\frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha}}{dt} = 0
+\frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha,\text{thermal}}}{dt} = 0
 $$
 
 Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions.
@@ -56,17 +56,6 @@ Where $\tau_p$ is the particle confinement time which we can assume is approxima
     The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1].
 
 
-
-
-### METIS Alpha Confinement
-
-$$
-\tau_{\alpha} = f_{\alpha}\tau_{\text{E}}\frac{R}{1-R}\tau_{\text{ne}}
-$$
-
-This is the model currently in METIS[^2] and is found in [^3]
-
-
 --------------
 
 ### Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()`
@@ -93,10 +82,10 @@ $$
 
 ---------------
 
-### Alpha Particle Flow Rate | `calculate_plasma_alphas_flow_rate()`
+### Thermal Alpha Particle Flow Rate | `calculate_plasma_alphas_thermal_flow_rate()`
 
 $$
-\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*}
+\frac{dn_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
 $$
 
 -----------------
@@ -177,11 +166,11 @@ This constraint ensures that the change in helium-3 particles as a function of t
 
 ----------------
 
-### Alpha Particle Flow Consistency
+### Thermal Alpha Particle Flow Consistency
 
 This constraint can be activated by stating `icc = 96` in the input file.
 
-This constraint ensures that the change in alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_flow_rate(()` is zero
+This constraint ensures that the change in thermal alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_thermal_flow_rate()` is zero
 
 **It is recommended to have this constraint on as it is a plasma consistency model**
 
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 7112e50496..65aa590d5f 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16632,7 +16632,7 @@ def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
         for j, f_t_alpha_energy_confinement in enumerate(
             f_t_alpha_energy_confinement_range
         ):
-            alpha_flow[i, j] = PlasmaFuelling.calculate_plasma_alphas_flow_rate(
+            alpha_flow[i, j] = PlasmaFuelling.calculate_plasma_alphas_thermal_flow_rate(
                 fusrat_dt_total=fusion_dt,
                 fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
                 t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index d6be3e39e5..47a7a5c6ef 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2059,18 +2059,18 @@ def constraint_equation_95(constraint_registration, data):
 @ConstraintManager.register_constraint(96, "particles/s", "=")
 def constraint_equation_96(constraint_registration, data):
     """
-    Alpha particle balance consistency equation.
+    Thermal alpha particle balance consistency equation.
 
-    This constraint ensures that the alpha particle balance is maintained in the
+    This constraint ensures that the thermal alpha particle balance is maintained in the
     plasma.
 
-    The numerator represents the total alpha particle sources,
-    while the denominator represents the total alpha particle sinks, which are
+    The numerator represents the total thermal alpha particle sources,
+    while the denominator represents the total thermal alpha particle sinks, which are
     related to the plasma volume, fuel ion density, and confinement time.
 
-    The equation ensures that the rate of alpha particle production equals the rate of
-    alpha particle loss, maintaining a steady-state condition for alpha particles in
-    the plasma.
+    The equation ensures that the rate of thermal alpha particle production equals the
+    rate of thermal alpha particle loss, maintaining a steady-state condition for
+    thermal alpha particles in the plasma.
     """
     # Alpha particle balance
     numerator = data.physics.fusrat_dt_total + data.physics.fusrat_plasma_dhe3
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 111df2e745..6ad8af5b39 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -555,7 +555,7 @@ def calculate_plasma_helium3_flow_rate(
         )
 
     @staticmethod
-    def calculate_plasma_alphas_flow_rate(
+    def calculate_plasma_alphas_thermal_flow_rate(
         fusrat_dt_total: float,
         fusrat_plasma_dhe3: float,
         t_energy_confinement: float,
@@ -563,7 +563,7 @@ def calculate_plasma_alphas_flow_rate(
         nd_plasma_alphas_thermal_vol_avg: float,
         vol_plasma: float,
     ) -> float:
-        """Calculate the net alpha particle flow rate into the plasma.
+        """Calculate the net thermal alpha particle flow rate into the plasma.
 
         Parameters
         ----------
@@ -583,16 +583,15 @@ def calculate_plasma_alphas_flow_rate(
         Returns
         -------
         float
-            Alpha particle flow rate into the plasma (particles/s).
+            Thermal alpha particle flow rate into the plasma (particles/s).
 
         Notes
         -----
-        - A positive value indicates a net flow of alpha particles into the plasma,
-        while a negative value indicates a net loss of alpha particles from the plasma.
+        - A positive value indicates a net flow of thermal alpha particles into the
+        plasma, while a negative value indicates a net loss of thermal alpha particles
+        from the plasma.
 
         """
-        # Alpha particle balance
-
         return (
             fusrat_dt_total
             + fusrat_plasma_dhe3

From a760073526a26dd2c2780f187ae1077e877546c8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 09:53:27 +0100
Subject: [PATCH 25/39] Refactor plasma fuelling calculations to include
 separate methods for thermal alpha particle source and loss rates; update
 related calculations in constraints and plotting functions.

---
 process/core/io/plot/summary.py    | 20 +++++-----
 process/core/solver/constraints.py | 17 +++++---
 process/models/physics/fuelling.py | 64 +++++++++++++++++++++++++++---
 3 files changed, 81 insertions(+), 20 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 65aa590d5f..5aefb6cc2d 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16632,15 +16632,17 @@ def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
         for j, f_t_alpha_energy_confinement in enumerate(
             f_t_alpha_energy_confinement_range
         ):
-            alpha_flow[i, j] = PlasmaFuelling.calculate_plasma_alphas_thermal_flow_rate(
-                fusrat_dt_total=fusion_dt,
-                fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
-                t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
-                nd_plasma_alphas_thermal_vol_avg=mfile.get(
-                    "nd_plasma_alphas_thermal_vol_avg", scan=scan
-                ),
-                vol_plasma=mfile.get("vol_plasma", scan=scan),
-                f_t_alpha_energy_confinement=f_t_alpha_energy_confinement,
+            alpha_flow[i, j] = (
+                PlasmaFuelling().calculate_plasma_alphas_thermal_flow_rate(
+                    fusrat_dt_total=fusion_dt,
+                    fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
+                    t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
+                    nd_plasma_alphas_thermal_vol_avg=mfile.get(
+                        "nd_plasma_alphas_thermal_vol_avg", scan=scan
+                    ),
+                    vol_plasma=mfile.get("vol_plasma", scan=scan),
+                    f_t_alpha_energy_confinement=f_t_alpha_energy_confinement,
+                )
             )
 
     contour = axis.contourf(
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 47a7a5c6ef..cf49285942 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2072,13 +2072,18 @@ def constraint_equation_96(constraint_registration, data):
     rate of thermal alpha particle loss, maintaining a steady-state condition for
     thermal alpha particles in the plasma.
     """
-    # Alpha particle balance
-    numerator = data.physics.fusrat_dt_total + data.physics.fusrat_plasma_dhe3
-    denominator = (
-        data.physics.nd_plasma_alphas_thermal_vol_avg * data.physics.vol_plasma
-    ) / (data.physics.t_energy_confinement * data.physics.f_t_alpha_energy_confinement)
+    source = PlasmaFuelling.calculate_plasma_alphas_thermal_source_rate(
+        fusrat_dt_total=data.physics.fusrat_dt_total,
+        fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3,
+    )
+    sink = PlasmaFuelling.calculate_plasma_alphas_thermal_loss_rate(
+        nd_plasma_alphas_thermal_vol_avg=data.physics.nd_plasma_alphas_thermal_vol_avg,
+        vol_plasma=data.physics.vol_plasma,
+        t_energy_confinement=data.physics.t_energy_confinement,
+        f_t_alpha_energy_confinement=data.physics.f_t_alpha_energy_confinement,
+    )
 
-    return eq(numerator, denominator, constraint_registration)
+    return eq(source, -sink, constraint_registration)
 
 
 @ConstraintManager.register_constraint(97, "", "=")
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 6ad8af5b39..1ee76a805d 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -555,7 +555,59 @@ def calculate_plasma_helium3_flow_rate(
         )
 
     @staticmethod
+    def calculate_plasma_alphas_thermal_source_rate(
+        fusrat_dt_total: float,
+        fusrat_plasma_dhe3: float,
+    ) -> float:
+        """Calculate the thermal alpha particle source rate in the plasma.
+
+        Parameters
+        ----------
+        fusrat_dt_total : float
+            Total DT fusion rate (particles/s).
+        fusrat_plasma_dhe3 : float
+            Deuterium consumption rate from D-He3 fusion (particles/s).
+
+        Returns
+        -------
+        float
+            Thermal alpha particle source rate in the plasma (particles/s).
+
+        """
+        return fusrat_dt_total + fusrat_plasma_dhe3
+
+    @staticmethod
+    def calculate_plasma_alphas_thermal_loss_rate(
+        t_energy_confinement: float,
+        f_t_alpha_energy_confinement: float,
+        nd_plasma_alphas_thermal_vol_avg: float,
+        vol_plasma: float,
+    ) -> float:
+        """Calculate the thermal alpha particle loss rate from the plasma.
+
+        Parameters
+        ----------
+        t_energy_confinement : float
+            Energy confinement time (s).
+        f_t_alpha_energy_confinement : float
+            Ratio of alpha particle confinement time to energy confinement time (dimensionless).
+        nd_plasma_alphas_thermal_vol_avg : float
+            Volume-averaged density of thermal alpha particles in the plasma (particles/m³).
+        vol_plasma : float
+            Plasma volume (m³).
+
+        Returns
+        -------
+        float
+            Thermal alpha particle loss rate from the plasma (particles/s).
+
+        """
+        return -(nd_plasma_alphas_thermal_vol_avg * vol_plasma) / (
+            t_energy_confinement * f_t_alpha_energy_confinement
+        )
+
     def calculate_plasma_alphas_thermal_flow_rate(
+        self,
         fusrat_dt_total: float,
         fusrat_plasma_dhe3: float,
         t_energy_confinement: float,
@@ -592,11 +644,13 @@ def calculate_plasma_alphas_thermal_flow_rate(
         from the plasma.
 
         """
-        return (
-            fusrat_dt_total
-            + fusrat_plasma_dhe3
-            - (nd_plasma_alphas_thermal_vol_avg * vol_plasma)
-            / (t_energy_confinement * f_t_alpha_energy_confinement)
+        return self.calculate_plasma_alphas_thermal_source_rate(
+            fusrat_dt_total=fusrat_dt_total, fusrat_plasma_dhe3=fusrat_plasma_dhe3
+        ) + self.calculate_plasma_alphas_thermal_loss_rate(
+            t_energy_confinement=t_energy_confinement,
+            f_t_alpha_energy_confinement=f_t_alpha_energy_confinement,
+            nd_plasma_alphas_thermal_vol_avg=nd_plasma_alphas_thermal_vol_avg,
+            vol_plasma=vol_plasma,
         )
 
     def output_fuelling_info(self):

From a99cc444fbc3ea3cb13840a258ca53147e5c27bc Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 10:26:32 +0100
Subject: [PATCH 26/39] Add beam deuterium and tritium injection rate variables
 to PhysicsData class

---
 process/data_structure/physics_variables.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index a51ad2325e..34251a3946 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1377,6 +1377,12 @@ class PhysicsData:
     f_molflow_plasma_fuelling_helium3: float = 0.0
     """Fraction of plasma fuelling that is helium-3"""
 
+    molflow_beam_deuterium_vv_injected: float = 0.0
+    """Beam deuterium injection rate into the vacuum vessel [particles/s]"""
+
+    molflow_beam_tritium_vv_injected: float = 0.0
+    """Beam tritium injection rate into the vacuum vessel [particles/s]"""
+
     q95_min: float = 0.0
     """Plasmalower limit for edge safety factor"""
 

From cd0e858898236d9c39b0e414b7cd5ddb7ad3a6e3 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 11:05:37 +0100
Subject: [PATCH 27/39] Add beam deuterium and tritium injection rate
 calculations to PlasmaFuelling and update related constraints and plotting
 functions

---
 .../source/physics-models/plasma_fuelling.md  | 20 ++---
 process/core/io/plot/summary.py               | 18 +++-
 process/core/solver/constraints.py            |  2 +
 process/models/physics/fuelling.py            | 88 ++++++++++++++++---
 process/models/physics/physics.py             | 11 +++
 5 files changed, 111 insertions(+), 28 deletions(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index 3a49869334..15046c7fa7 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -5,11 +5,11 @@
 The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the thermal helium $\alpha$ ash.
 
 $$
-\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
 $$
 
 $$
-\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
 $$
 
 $$
@@ -26,7 +26,7 @@ $$
 \frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha,\text{thermal}}}{dt} = 0
 $$
 
-Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions.
+Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions. $\Gamma_{\text{T,beam}}$ and $\Gamma_{\text{D,beam}}$ respectively are the particle source rates from neutral beam systems (if present).
 
 
 The fuelling fractional compositions is given by $f$
@@ -61,7 +61,7 @@ Where $\tau_p$ is the particle confinement time which we can assume is approxima
 ### Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()`
 
 $$
-\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
 $$
 
 ---------------
@@ -69,7 +69,7 @@ $$
 ### Deuterium Flow Rate | `calculate_plasma_deuterium_flow_rate()`
 
 $$
-\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
 $$
 
 ---------------
@@ -95,10 +95,10 @@ $$
 In a steady state tokamak, the burnup fraction ($f_b$) is explicitly defined by the following rate equation:
 
 $$
-f_b = \frac{\Gamma_{\text{fusion}}}{\Gamma_{\text{fuel}}}
+f_b = \frac{\Gamma_{\text{fusion}}}{\sum\Gamma_{\text{fuel}}}
 $$
 
-where $\Gamma_{\text{fusion}}$ is the fusion reaction rate per second $[\text{s}^{-1}]$, and $\Gamma_{\text{fuel}}$ is the fuel injection rate into the vacuum vessel in particles per second $[\text{s}^{-1}]$.
+where $\Gamma_{\text{fusion}}$ is the fusion reaction rate per second $[\text{s}^{-1}]$, and $\sum\Gamma_{\text{fuel}}$ is the sum of all fuel injection rates into the vacuum vessel by any means in particles per second $[\text{s}^{-1}]$.
 
 -----------------
 
@@ -107,7 +107,7 @@ where $\Gamma_{\text{fusion}}$ is the fusion reaction rate per second $[\text{s}
 For the total burnup fraction we state:
 
 $$
-\overbrace{f_b}^{\texttt{f_plasma_fuel_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+T}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}}
+\overbrace{f_b}^{\texttt{f_plasma_fuel_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+T}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}}+\Gamma_{\text{T,beam}}}
 $$
 
 Here the factor of 2 is included as each fusion reaction removes 2 particles but our fuelling rate looks at indivudal particles injected.
@@ -119,7 +119,7 @@ Here the factor of 2 is included as each fusion reaction removes 2 particles but
 For just the tritium burnup fraction we state:
 
 $$
-\overbrace{f_b}^{\texttt{f_plasma_tritium_burnup}} = \frac{\left(\Gamma_{\text{D+T}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,T}}} 
+\overbrace{f_b}^{\texttt{f_plasma_tritium_burnup}} = \frac{\left(\Gamma_{\text{D+T}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,T}}+\Gamma_{\text{T,beam}}} 
 $$
 
 ------------------
@@ -129,7 +129,7 @@ $$
 For just the deuterium burnup fraction we state:
 
 $$
-\overbrace{f_b}^{\texttt{f_plasma_deuterium_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,D}}}
+\overbrace{f_b}^{\texttt{f_plasma_deuterium_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,D}}+\Gamma_{\text{D,beam}}}
 $$
 
 ------------------
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 5aefb6cc2d..ee9c1e50ce 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16500,6 +16500,9 @@ def plot_tritium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
                 molflow_plasma_fuelling_vv_injected=mfile.get(
                     "molflow_plasma_fuelling_vv_injected", scan=scan
                 ),
+                molflow_beam_tritium_vv_injected=mfile.get(
+                    "molflow_beam_tritium_vv_injected", scan=scan
+                ),
                 fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan),
                 fusrat_plasma_dd_triton=mfile.get("fusrat_plasma_dd_triton", scan=scan),
                 t_energy_confinement=mfile.get("t_energy_confinement", scan=scan),
@@ -16568,6 +16571,9 @@ def plot_deuterium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int):
                 molflow_plasma_fuelling_vv_injected=mfile.get(
                     "molflow_plasma_fuelling_vv_injected", scan=scan
                 ),
+                molflow_beam_deuterium_vv_injected=mfile.get(
+                    "molflow_beam_deuterium_vv_injected", scan=scan
+                ),
                 fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan),
                 fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan),
                 fusrat_plasma_dd_total=mfile.get("fusrat_plasma_dd_total", scan=scan),
@@ -16755,13 +16761,17 @@ def plot_fuelling_info(fig: plt.Figure, mfile: MFile, scan: int):
     """Plot fuelling information."""
     msg = (
         f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n"
-        f"Total fuelling rate:"
+        f"Fuelling rate into vessel:"
         f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n"
-        f"Total fuelling rate: "
+        f"Fuelling rate into vessel: "
         f"{mfile.get('molflow_plasma_fuelling_vv_injected_moles', scan=scan):.4e} moles/s\n"
-        f"Total fuelling loss: "
+        f"Additional deuterium beam fuelling: "
+        f"{mfile.get('molflow_beam_deuterium_vv_injected', scan=scan):.4e} particles/s\n"
+        f"Additional tritium beam fuelling: "
+        f"{mfile.get('molflow_beam_tritium_vv_injected', scan=scan):.4e} particles/s\n\n"
+        f"Fuelling loss that doesn't enter the plasma: "
         f"{mfile.get('molflow_plasma_fuelling_loss', scan=scan):.4e} particles/s\n"
-        f"Total fuelling loss: "
+        f"Fuelling loss that doesn't enter the plasma: "
         f"{mfile.get('molflow_plasma_fuelling_loss_moles', scan=scan):.4e} moles/s\n"
         f"Fuelling Rate Efficiency ($\\eta_{{\\text{{fuelling}}}}$): "
         f"{mfile.get('eta_plasma_fuelling', scan=scan):.4f}\n"
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index cf49285942..74c9dccfe7 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -1976,6 +1976,7 @@ def constraint_equation_93(constraint_registration, data):
         eta_plasma_fuelling=data.physics.eta_plasma_fuelling,
         molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected,
         fusrat_plasma_dd_triton=data.physics.fusrat_plasma_dd_triton,
+        molflow_beam_tritium_vv_injected=data.physics.molflow_beam_tritium_vv_injected,
     )
     sink = PlasmaFuelling.calculate_plasma_tritium_loss_rate(
         fusrat_dt_total=data.physics.fusrat_dt_total,
@@ -2008,6 +2009,7 @@ def constraint_equation_94(constraint_registration, data):
         f_molflow_plasma_fuelling_deuterium=data.physics.f_molflow_plasma_fuelling_deuterium,
         eta_plasma_fuelling=data.physics.eta_plasma_fuelling,
         molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected,
+        molflow_beam_deuterium_vv_injected=data.physics.molflow_beam_deuterium_vv_injected,
     )
     sink = PlasmaFuelling.calculate_plasma_deuterium_loss_rate(
         fusrat_dt_total=data.physics.fusrat_dt_total,
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 1ee76a805d..6dd1f776d7 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -31,16 +31,20 @@ def run(self):
         self.data.physics.f_plasma_fuel_burnup = self.calculate_fuel_burnup_fraction(
             fusrat_total=self.data.physics.fusrat_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+            molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
+            molflow_beam_tritium_vv_injected=self.data.physics.molflow_beam_tritium_vv_injected,
         )
         self.data.physics.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction(
             fusrat_dt_total=self.data.physics.fusrat_dt_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
             f_molflow_plasma_fuelling_tritium=self.data.physics.f_molflow_plasma_fuelling_tritium,
+            molflow_beam_tritium_vv_injected=self.data.physics.molflow_beam_tritium_vv_injected,
         )
 
         self.data.physics.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction(
             fusrat_plasma_dd_total=self.data.physics.fusrat_plasma_dd_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+            molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
             f_molflow_plasma_fuelling_deuterium=self.data.physics.f_molflow_plasma_fuelling_deuterium,
             fusrat_dt_total=self.data.physics.fusrat_dt_total,
             fusrat_plasma_dhe3=self.data.physics.fusrat_plasma_dhe3,
@@ -53,7 +57,10 @@ def output(self):
 
     @staticmethod
     def calculate_fuel_burnup_fraction(
-        fusrat_total: float, molflow_plasma_fuelling_vv_injected: float
+        fusrat_total: float,
+        molflow_plasma_fuelling_vv_injected: float,
+        molflow_beam_deuterium_vv_injected: float,
+        molflow_beam_tritium_vv_injected: float,
     ) -> float:
         """Calculate the fuel burnup fraction
 
@@ -63,6 +70,10 @@ def calculate_fuel_burnup_fraction(
             Total fusion rate (particles/s).
         molflow_plasma_fuelling_vv_injected : float
             Total fuelling rate into vacuum vessel (particles/s).
+        molflow_beam_deuterium_vv_injected : float
+            Beam deuterium injection rate into the vacuum vessel (particles/s).
+        molflow_beam_tritium_vv_injected : float
+            Beam tritium injection rate into the vacuum vessel (particles/s).
 
         Returns
         -------
@@ -74,13 +85,22 @@ def calculate_fuel_burnup_fraction(
         as the fuelling rate is in particles/s.
 
         """
-        return 2 * fusrat_total / molflow_plasma_fuelling_vv_injected
+        return (
+            2
+            * fusrat_total
+            / (
+                molflow_plasma_fuelling_vv_injected
+                + molflow_beam_deuterium_vv_injected
+                + molflow_beam_tritium_vv_injected
+            )
+        )
 
     @staticmethod
     def calculate_tritium_burnup_fraction(
         fusrat_dt_total: float,
         molflow_plasma_fuelling_vv_injected: float,
         f_molflow_plasma_fuelling_tritium: float,
+        molflow_beam_tritium_vv_injected: float,
     ) -> float:
         """Calculate the tritium burnup fraction
 
@@ -92,6 +112,8 @@ def calculate_tritium_burnup_fraction(
             Total fuelling rate into vacuum vessel (particles/s).
         f_molflow_plasma_fuelling_tritium : float
             Fraction of tritium in the plasma fuelling.
+        molflow_beam_tritium_vv_injected : float
+            Beam tritium injection rate into the vacuum vessel (particles/s).
 
         Returns
         -------
@@ -104,13 +126,15 @@ def calculate_tritium_burnup_fraction(
 
         """
         return fusrat_dt_total / (
-            molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_tritium
+            (molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_tritium)
+            + molflow_beam_tritium_vv_injected
         )
 
     @staticmethod
     def calculate_deuterium_burnup_fraction(
         fusrat_dt_total: float,
         molflow_plasma_fuelling_vv_injected: float,
+        molflow_beam_deuterium_vv_injected: float,
         f_molflow_plasma_fuelling_deuterium: float,
         fusrat_plasma_dd_total: float,
         fusrat_plasma_dhe3: float,
@@ -123,6 +147,8 @@ def calculate_deuterium_burnup_fraction(
             Total DT fusion rate (particles/s).
         molflow_plasma_fuelling_vv_injected : float
             Total fuelling rate into vacuum vessel (particles/s).
+        molflow_beam_deuterium_vv_injected : float
+            Beam deuterium injection rate into the vacuum vessel (particles/s).
         f_molflow_plasma_fuelling_deuterium : float
             Fraction of deuterium in the plasma fuelling.
         fusrat_plasma_dd_total : float
@@ -141,7 +167,8 @@ def calculate_deuterium_burnup_fraction(
 
         """
         return (fusrat_dt_total + 2 * fusrat_plasma_dd_total + fusrat_plasma_dhe3) / (
-            molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium
+            (molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium)
+            + molflow_beam_deuterium_vv_injected
         )
 
     def calculate_plasma_tritium_flow_rate(
@@ -149,6 +176,7 @@ def calculate_plasma_tritium_flow_rate(
         f_molflow_plasma_fuelling_tritium: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
+        molflow_beam_tritium_vv_injected: float,
         fusrat_dt_total: float,
         fusrat_plasma_dd_triton: float,
         t_energy_confinement: float,
@@ -167,6 +195,8 @@ def calculate_plasma_tritium_flow_rate(
             Fuelling rate efficiency.
         molflow_plasma_fuelling_vv_injected : float
             Total fuelling rate (particles/s).
+        molflow_beam_tritium_vv_injected : float
+            Beam tritium injection rate into the vacuum vessel (particles/s).
         fusrat_dt_total : float
             Total DT fusion rate (particles/s).
         fusrat_plasma_dd_triton : float
@@ -198,6 +228,7 @@ def calculate_plasma_tritium_flow_rate(
             eta_plasma_fuelling=eta_plasma_fuelling,
             molflow_plasma_fuelling_vv_injected=molflow_plasma_fuelling_vv_injected,
             fusrat_plasma_dd_triton=fusrat_plasma_dd_triton,
+            molflow_beam_tritium_vv_injected=molflow_beam_tritium_vv_injected,
         ) + self.calculate_plasma_tritium_loss_rate(
             fusrat_dt_total=fusrat_dt_total,
             t_energy_confinement=t_energy_confinement,
@@ -213,6 +244,7 @@ def calculate_plasma_tritium_source_rate(
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
         fusrat_plasma_dd_triton: float,
+        molflow_beam_tritium_vv_injected: float,
     ) -> float:
         """Calculate the tritium source rate in the plasma.
 
@@ -226,6 +258,8 @@ def calculate_plasma_tritium_source_rate(
             Total fuelling rate (particles/s).
         fusrat_plasma_dd_triton : float
             Tritium production rate from D-D fusion (particles/s).
+        molflow_beam_tritium_vv_injected : float
+            Beam tritium injection rate into the vacuum vessel (particles/s).
 
         Returns
         -------
@@ -234,10 +268,14 @@ def calculate_plasma_tritium_source_rate(
 
         """
         return (
-            f_molflow_plasma_fuelling_tritium
-            * eta_plasma_fuelling
-            * molflow_plasma_fuelling_vv_injected
-        ) + fusrat_plasma_dd_triton
+            (
+                f_molflow_plasma_fuelling_tritium
+                * eta_plasma_fuelling
+                * molflow_plasma_fuelling_vv_injected
+            )
+            + fusrat_plasma_dd_triton
+            + molflow_beam_tritium_vv_injected
+        )
 
     @staticmethod
     def calculate_plasma_tritium_loss_rate(
@@ -281,6 +319,7 @@ def calculate_plasma_deuterium_source_rate(
         f_molflow_plasma_fuelling_deuterium: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
+        molflow_beam_deuterium_vv_injected: float,
     ) -> float:
         """Calculate the deuterium source rate in the plasma.
 
@@ -292,6 +331,8 @@ def calculate_plasma_deuterium_source_rate(
             Fuelling rate efficiency.
         molflow_plasma_fuelling_vv_injected : float
             Total fuelling rate (particles/s).
+        molflow_beam_deuterium_vv_injected : float
+            Beam deuterium injection rate into the vacuum vessel (particles/s).
 
         Returns
         -------
@@ -303,7 +344,7 @@ def calculate_plasma_deuterium_source_rate(
             f_molflow_plasma_fuelling_deuterium
             * eta_plasma_fuelling
             * molflow_plasma_fuelling_vv_injected
-        )
+        ) + molflow_beam_deuterium_vv_injected
 
     @staticmethod
     def calculate_plasma_deuterium_loss_rate(
@@ -358,6 +399,7 @@ def calculate_plasma_deuterium_flow_rate(
         f_molflow_plasma_fuelling_deuterium: float,
         eta_plasma_fuelling: float,
         molflow_plasma_fuelling_vv_injected: float,
+        molflow_beam_deuterium_vv_injected: float,
         fusrat_dt_total: float,
         fusrat_plasma_dhe3: float,
         fusrat_plasma_dd_total: float,
@@ -377,6 +419,8 @@ def calculate_plasma_deuterium_flow_rate(
             Fuelling rate efficiency.
         molflow_plasma_fuelling_vv_injected : float
             Total fuelling rate (particles/s).
+        molflow_beam_deuterium_vv_injected : float
+            Beam deuterium injection rate into the vacuum vessel (particles/s).
         fusrat_dt_total : float
             Total DT fusion rate (particles/s).
         fusrat_plasma_dhe3 : float
@@ -410,6 +454,7 @@ def calculate_plasma_deuterium_flow_rate(
             f_molflow_plasma_fuelling_deuterium=f_molflow_plasma_fuelling_deuterium,
             eta_plasma_fuelling=eta_plasma_fuelling,
             molflow_plasma_fuelling_vv_injected=molflow_plasma_fuelling_vv_injected,
+            molflow_beam_deuterium_vv_injected=molflow_beam_deuterium_vv_injected,
         ) + self.calculate_plasma_deuterium_loss_rate(
             fusrat_dt_total=fusrat_dt_total,
             fusrat_plasma_dd_total=fusrat_plasma_dd_total,
@@ -658,33 +703,48 @@ def output_fuelling_info(self):
         po.oheadr(self.outfile, "Plasma Fuelling")
         po.ovarre(
             self.outfile,
-            "Fuelling rate (nucleus-pairs/s)",
+            "Fuelling rate into vessel [particles/s]",
             "(molflow_plasma_fuelling_vv_injected)",
             self.data.physics.molflow_plasma_fuelling_vv_injected,
             "OP ",
         )
         po.ovarre(
             self.outfile,
-            "Fuelling rate (moles/s)",
+            "Fuelling rate into vessel [moles/s]",
             "(molflow_plasma_fuelling_vv_injected_moles)",
             self.data.physics.molflow_plasma_fuelling_vv_injected_moles,
             "OP ",
         )
         po.ovarre(
             self.outfile,
-            "Fuelling loss (nucleus-pairs/s)",
+            "Fuelling loss that doesn't enter the plasma [particles/s]",
             "(molflow_plasma_fuelling_loss)",
             self.data.physics.molflow_plasma_fuelling_loss,
             "OP ",
         )
         po.ovarre(
             self.outfile,
-            "Fuelling loss (moles/s)",
+            "Fuelling loss that doesn't enter the plasma [moles/s]",
             "(molflow_plasma_fuelling_loss_moles)",
             self.data.physics.molflow_plasma_fuelling_loss_moles,
             "OP ",
         )
         po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Beam deuterium injection rate into the vacuum vessel [particles/s]",
+            "(molflow_beam_deuterium_vv_injected)",
+            self.data.physics.molflow_beam_deuterium_vv_injected,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Beam tritium injection rate into the vacuum vessel [particles/s]",
+            "(molflow_beam_tritium_vv_injected)",
+            self.data.physics.molflow_beam_tritium_vv_injected,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "Fraction of plasma fuelling that is deuterium",
@@ -723,7 +783,7 @@ def output_fuelling_info(self):
         )
         po.ovarre(
             self.outfile,
-            "Fuel burn-up rate (reactions/s)",
+            "Fuel burn-up rate [reactions/s]",
             "(fusrat_total)",
             self.data.physics.fusrat_total,
             "OP ",
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 79ef1f58cc..7c5aa37512 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -672,6 +672,17 @@ def run(self):
             self.data.physics.p_beam_dt_mw = self.data.physics.p_beam_alpha_mw * (
                 1 / (1 - constants.DT_NEUTRON_ENERGY_FRACTION)
             )
+
+            self.data.physics.molflow_beam_deuterium_vv_injected = (
+                self.data.current_drive.c_beam_total
+                * (1.0 - self.data.current_drive.f_beam_tritium)
+                / constants.ELECTRON_CHARGE
+            )
+            self.data.physics.molflow_beam_tritium_vv_injected = (
+                self.data.current_drive.c_beam_total
+                * self.data.current_drive.f_beam_tritium
+                / constants.ELECTRON_CHARGE
+            )
         else:
             # If no beams present then the total alpha rates and power are the same as
             # the plasma values

From 60b523fbbda9232ff8506e38767c232e847cd8d0 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 13:15:48 +0100
Subject: [PATCH 28/39] Fix particle balance equations in Plasma Fuelling
 documentation to use consistent notation for time derivatives.

---
 .../source/physics-models/plasma_fuelling.md   | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index 15046c7fa7..00124de35d 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -5,25 +5,25 @@
 The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the thermal helium $\alpha$ ash.
 
 $$
-\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+\frac{dN_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
 $$
 
 $$
-\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+\frac{dN_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
 $$
 
 $$
-\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} -\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+\frac{dN_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} -\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
 $$
 
 $$
-\frac{dn_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
+\frac{dN_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
 $$
 
 In a steady state equilibrium all 4 of these equations should balance, therefore:
 
 $$
-\frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha,\text{thermal}}}{dt} = 0
+\frac{dN_{\text{D}}}{dt} = \frac{dN_{\text{T}}}{dt} = \frac{dN_{\text{3He}}}{dt} = \frac{dN_{\alpha,\text{thermal}}}{dt} = 0
 $$
 
 Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions. $\Gamma_{\text{T,beam}}$ and $\Gamma_{\text{D,beam}}$ respectively are the particle source rates from neutral beam systems (if present).
@@ -61,7 +61,7 @@ Where $\tau_p$ is the particle confinement time which we can assume is approxima
 ### Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()`
 
 $$
-\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+\frac{dN_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
 $$
 
 ---------------
@@ -69,7 +69,7 @@ $$
 ### Deuterium Flow Rate | `calculate_plasma_deuterium_flow_rate()`
 
 $$
-\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+\frac{dN_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
 $$
 
 ---------------
@@ -77,7 +77,7 @@ $$
 ### Helium-3 Flow Rate | `calculate_plasma_helium3_flow_rate()`
 
 $$
-\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}}-\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+\frac{dN_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}}-\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
 $$
 
 ---------------
@@ -85,7 +85,7 @@ $$
 ### Thermal Alpha Particle Flow Rate | `calculate_plasma_alphas_thermal_flow_rate()`
 
 $$
-\frac{dn_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
+\frac{dN_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
 $$
 
 -----------------

From 432a9f80c75340869e43d57b8e08fd0c7cc02dfe Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 23 Jul 2026 13:35:58 +0100
Subject: [PATCH 29/39] Update plasma fuelling documentation to clarify
 constraints as plasma equilibrium solution models

---
 documentation/source/physics-models/plasma_fuelling.md | 10 ++++++----
 process/core/solver/constraints.py                     |  4 ++--
 process/models/physics/fuelling.py                     |  9 ++++++---
 process/models/stellarator/stellarator.py              |  2 ++
 4 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index 00124de35d..488656faf5 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -136,13 +136,15 @@ $$
 
 ## Key Constraints
 
+The implementation of the equations above as equality constraints is required in order to achieve a plasma steady state equilibrium where there is not a net rate of change of species number.
+
 ### Tritium Flow Consistency
 
 This constraint can be activated by stating `icc = 93` in the input file.
 
 This constraint ensures that the change in tritium particles as a function of time is zero. It ensures the output of `calculate_plasma_tritium_flow_rate()` is zero
 
-**It is recommended to have this constraint on as it is a plasma consistency model**
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
 
 -----------------
 
@@ -152,7 +154,7 @@ This constraint can be activated by stating `icc = 94` in the input file.
 
 This constraint ensures that the change in deuterium particles as a function of time is zero. It ensures the output of `calculate_plasma_deuterium_flow_rate()` is zero
 
-**It is recommended to have this constraint on as it is a plasma consistency model**
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
 
 ----------------
 
@@ -162,7 +164,7 @@ This constraint can be activated by stating `icc = 95` in the input file.
 
 This constraint ensures that the change in helium-3 particles as a function of time is zero. It ensures the output of `calculate_plasma_helium3_flow_rate()` is zero
 
-**It is recommended to have this constraint on as it is a plasma consistency model**
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
 
 ----------------
 
@@ -172,7 +174,7 @@ This constraint can be activated by stating `icc = 96` in the input file.
 
 This constraint ensures that the change in thermal alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_thermal_flow_rate()` is zero
 
-**It is recommended to have this constraint on as it is a plasma consistency model**
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
 
 ------------------
 
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 74c9dccfe7..6826204d2c 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -2002,8 +2002,8 @@ def constraint_equation_94(constraint_registration, data):
     reactions, while the denominator represents the total deuterium sinks, which are
     related to the plasma volume, fuel ion density, and confinement time.
 
-    The equation ensures that the rate of deuterium production equals the rate of deuterium
-    loss, maintaining a steady-state condition for deuterium in the plasma.
+    The equation ensures that the rate of deuterium production equals the rate of
+    deuterium loss, maintaining a steady-state condition for deuterium in the plasma.
     """
     source = PlasmaFuelling.calculate_plasma_deuterium_source_rate(
         f_molflow_plasma_fuelling_deuterium=data.physics.f_molflow_plasma_fuelling_deuterium,
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 6dd1f776d7..3df85b9e99 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -635,9 +635,11 @@ def calculate_plasma_alphas_thermal_loss_rate(
         t_energy_confinement : float
             Energy confinement time (s).
         f_t_alpha_energy_confinement : float
-            Ratio of alpha particle confinement time to energy confinement time (dimensionless).
+            Ratio of alpha particle confinement time to energy confinement time
+            (dimensionless).
         nd_plasma_alphas_thermal_vol_avg : float
-            Volume-averaged density of thermal alpha particles in the plasma (particles/m³).
+            Volume-averaged density of thermal alpha particles in the plasma
+            (particles/m³).
         vol_plasma : float
             Plasma volume (m³).
 
@@ -671,7 +673,8 @@ def calculate_plasma_alphas_thermal_flow_rate(
         t_energy_confinement : float
             Energy confinement time (s).
         f_t_alpha_energy_confinement : float
-            Ratio of alpha particle confinement time to energy confinement time (dimensionless).
+            Ratio of alpha particle confinement time to energy confinement time
+            (dimensionless).
         nd_plasma_alphas_thermal_vol_avg : float
             Volume-averaged density of alpha particles in the plasma (particles/m³).
         vol_plasma : float
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 05d004deb2..3594636415 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2388,6 +2388,8 @@ def st_phys(self, output):
         self.data.physics.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction(
             fusrat_total=self.data.physics.fusrat_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
+            molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
+            molflow_beam_tritium_vv_injected=self.data.physics.molflow_beam_tritium_vv_injected,
         )
 
         # Calculate the neoclassical sanity check with PROCESS parameters

From 24148cd0e35b3ea9541d30c848baf367f1773332 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 28 Jul 2026 16:14:36 +0100
Subject: [PATCH 30/39] Fix particle recycling and fuelling efficiency indices
 in large tokamak input file

---
 .../input_files/large_tokamak_nof.IN.DAT         | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/tests/regression/input_files/large_tokamak_nof.IN.DAT b/tests/regression/input_files/large_tokamak_nof.IN.DAT
index 7cfd5e5ced..4d61773ef7 100644
--- a/tests/regression/input_files/large_tokamak_nof.IN.DAT
+++ b/tests/regression/input_files/large_tokamak_nof.IN.DAT
@@ -191,27 +191,27 @@ t_burn_min = 7200.0
 ***********************
 
 * Particle recycling fraction
-ixc = 177
+ixc = 178
 f_plasma_particles_lcfs_recycled = 0.9
 
 * Plasma fuelling efficiecy
-ixc = 178
+ixc = 179
 eta_plasma_fuelling = 0.7
 
 * Injected VV fuelling rate
-ixc = 179
+ixc = 180
 molflow_plasma_fuelling_vv_injected = 5e21
-boundl(179) = 1e20
+boundl(180) = 1e20
 
 * Deuterium fuelling fraction
-ixc = 180
+ixc = 181
 f_molflow_plasma_fuelling_deuterium = 0.5
-boundl(180) = 0.4
+boundl(181) = 0.4
 
 * Tritium fuelling fraction
-ixc = 181
+ixc = 182
 f_molflow_plasma_fuelling_tritium = 0.5
-boundl(181) = 0.4
+boundl(182) = 0.4
 
 
 * b_plasma_toroidal_on_axis [T]

From 7e39dd1351796907a6a4c4a2e71267d17d572cd5 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 28 Jul 2026 16:25:06 +0100
Subject: [PATCH 31/39] Remove duplicate output for reaction rates

---
 process/models/physics/physics.py | 27 +++++----------------------
 1 file changed, 5 insertions(+), 22 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 7c5aa37512..62b3d4d21d 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -1816,6 +1816,7 @@ def outplas(self):
             self.data.physics.fusrat_total,
             "OP ",
         )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "D-T Fusion rate: total (reactions/sec)",
@@ -1830,7 +1831,7 @@ def outplas(self):
             self.data.physics.fusrat_plasma_dt,
             "OP ",
         )
-
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "D-D -> 3He Fusion rate: plasma (reactions/sec)",
@@ -1852,6 +1853,7 @@ def outplas(self):
             self.data.physics.fusrat_plasma_dd_total,
             "OP ",
         )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "D-3He Fusion rate: total (reactions/sec)",
@@ -1859,6 +1861,7 @@ def outplas(self):
             self.data.physics.fusrat_plasma_dhe3,
             "OP ",
         )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "Neutron production rate: total (particles/sec)",
@@ -1866,27 +1869,7 @@ def outplas(self):
             self.data.physics.fusrat_neutron_production_total,
             "OP ",
         )
-        po.ovarre(
-            self.outfile,
-            "D-3He Fusion rate: total (reactions/sec)",
-            "(fusrat_plasma_dhe3)",
-            self.data.physics.fusrat_plasma_dhe3,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Neutron production rate: total (particles/sec)",
-            "(fusrat_neutron_production_total)",
-            self.data.physics.fusrat_neutron_production_total,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "D-D Fusion rate: total (reactions/sec)",
-            "(fusrat_plasma_dd_total)",
-            self.data.physics.fusrat_plasma_dd_total,
-            "OP ",
-        )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "Fusion rate density: total (reactions/m³/sec)",

From c6846a38bf18a2122b63208372f8744deff779bb Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 28 Jul 2026 16:27:07 +0100
Subject: [PATCH 32/39] Fix neutron production total calculation in Physics
 model

---
 process/models/physics/physics.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 62b3d4d21d..16a4d0054b 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -702,14 +702,14 @@ def run(self):
             + self.data.physics.fusrat_plasma_dd_triton
         )
 
-        self.data.physics.fusrat_neutron_production_total = (
-            self.data.physics.fusrat_plasma_dd_helion + self.data.physics.fusrat_dt_total
-        )
-
         self.data.physics.fusrat_dt_total = (
             self.data.physics.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY)
         )
 
+        self.data.physics.fusrat_neutron_production_total = (
+            self.data.physics.fusrat_plasma_dd_helion + self.data.physics.fusrat_dt_total
+        )
+
         # Create some derived values and add beam contribution to fusion power
         (
             self.data.physics.pden_neutron_total_mw,

From 69b0039bf93254909563b0a90b01bbe221aeb630 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 30 Jul 2026 09:52:06 +0100
Subject: [PATCH 33/39] Add particle balance constraints and iteration vars to
 the st_regression test

---
 .../input_files/st_regression.IN.DAT          | 39 ++++++++++++++++++-
 1 file changed, 38 insertions(+), 1 deletion(-)

diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT
index 2a346db5ac..9f6e8aeb25 100644
--- a/tests/regression/input_files/st_regression.IN.DAT
+++ b/tests/regression/input_files/st_regression.IN.DAT
@@ -55,10 +55,12 @@ i_figure_merit = -5
 * DESCRIPTION:   Switch for Figure-of-Merit 
 * JUSTIFICATION: Selected optimisation
 
-n_equality_constraints = 3
+n_equality_constraints = 7
 * DESCRIPTION:   Number of equality constraints (the first n_equality_constraints icc = ... statements are considered equality)
 * JUSTIFICATION: 3 equality (consistency) equations
 
+
+
 *________________________________________________________________________*
 *------------------------------ *Key Variables* -------------------------*
 *‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾*
@@ -85,6 +87,18 @@ icc = 11
 * JUSTIFICATION: Consistency equations should always be on
 * VARIABLES:     rmajor, Rest calculated in-situ
 
+* Tritium particle balance
+icc = 93
+
+* Deuterium particle balance
+icc = 94
+
+* Alpha particle balance
+icc = 96
+
+* Fuelling composition consistency
+icc = 97
+
 *________________________________________________________________________*
 *------------------------------- *Plasma* -------------------------------*
 *‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾*
@@ -555,6 +569,29 @@ f_sync_reflect = 0.6
 * allowance for radiation. This is not recommended for power plant models.
 * JUSTIFICATION: Turned off, assuming default = 1
 
+* 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
+
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 icc = 17 

From 6445e01e89f2791ccc38f3530aaffb2d48563c56 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 30 Jul 2026 10:05:52 +0100
Subject: [PATCH 34/39] Enhance plasma fuelling documentation problem input
 example and optimisation explanation

---
 .../source/physics-models/plasma_fuelling.md  | 46 ++++++++++++++++++-
 process/models/physics/fuelling.py            |  2 +
 2 files changed, 47 insertions(+), 1 deletion(-)

diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md
index 488656faf5..9b2423ea28 100644
--- a/documentation/source/physics-models/plasma_fuelling.md
+++ b/documentation/source/physics-models/plasma_fuelling.md
@@ -136,7 +136,51 @@ $$
 
 ## Key Constraints
 
-The implementation of the equations above as equality constraints is required in order to achieve a plasma steady state equilibrium where there is not a net rate of change of species number.
+The implementation of the equations above as equality constraints is required in order to achieve a plasma steady state equilibrium where there is not a net rate of change of species number. In any given run several of these constraints and iterations variables will almost always be on.Below is the generic input used for a D-T only plasma where we have allowed all of the key solution variables to be iteration variables: 
+
+```python
+* 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
+```
+
+Note that the Helium-3 consistency equation (`icc= 95`) is not active as we are not actively fuelling and have no stated Helium-3 in the plasma composition. 
+
+If solving this as an evaluation problem it would be typical to fix the values of the reycling fraction, $R$ (`f_plasma_particles_lcfs_recycled`) and the fuelling efficiency $\eta_{\text{fuelling}}$ (`eta_plasma_fuelling`) as these represent higher order assumptions already about the plant design that `PROCESS` cannot model. This realistically leaves the total fuelling rate $\Gamma_{\text{fuel}}$ (`molflow_plasma_fuelling_vv_injected`) and the fuelling fractions $f_{\text{fuelling,D}}$,$f_{\text{fuelling,T}}$ (`f_molflow_plasma_fuelling_deuterium,f_molflow_plasma_fuelling_tritium`) to solve this consistency equations. This is ideal as the fuelling rates and fractions are ideal optimisation paramters as they can be directly controlled from the machine control room.
+
+------------
 
 ### Tritium Flow Consistency
 
diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 3df85b9e99..a58f18814c 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -1,3 +1,5 @@
+"""Plasma fuelling model."""
+
 from process.core import constants
 from process.core import process_output as po
 from process.core.model import Model

From d51682e5ae3b2265dcb0e3600aec2cf45d3673ff Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Fri, 31 Jul 2026 14:49:44 +0100
Subject: [PATCH 35/39] Apply new fuelling constraints to large tokamak eval
 regression test

---
 process/models/physics/fuelling.py            |  4 +--
 process/models/stellarator/stellarator.py     |  2 +-
 .../input_files/large_tokamak_eval.IN.DAT     | 33 ++++++++++++++++++-
 3 files changed, 35 insertions(+), 4 deletions(-)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index a58f18814c..3fe44fad05 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -36,14 +36,14 @@ def run(self):
             molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
             molflow_beam_tritium_vv_injected=self.data.physics.molflow_beam_tritium_vv_injected,
         )
-        self.data.physics.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction(
+        self.data.physics.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction(  # noqa: E501
             fusrat_dt_total=self.data.physics.fusrat_dt_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
             f_molflow_plasma_fuelling_tritium=self.data.physics.f_molflow_plasma_fuelling_tritium,
             molflow_beam_tritium_vv_injected=self.data.physics.molflow_beam_tritium_vv_injected,
         )
 
-        self.data.physics.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction(
+        self.data.physics.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction(  # noqa: E501
             fusrat_plasma_dd_total=self.data.physics.fusrat_plasma_dd_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
             molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 3594636415..aa56547567 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2385,7 +2385,7 @@ def st_phys(self, output):
             self.data.physics.t_energy_confinement,
         )
 
-        self.data.physics.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction(
+        self.data.physics.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction(  # noqa: E501
             fusrat_total=self.data.physics.fusrat_total,
             molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected,
             molflow_beam_deuterium_vv_injected=self.data.physics.molflow_beam_deuterium_vv_injected,
diff --git a/tests/regression/input_files/large_tokamak_eval.IN.DAT b/tests/regression/input_files/large_tokamak_eval.IN.DAT
index efc4c37dd2..f630942ff4 100644
--- a/tests/regression/input_files/large_tokamak_eval.IN.DAT
+++ b/tests/regression/input_files/large_tokamak_eval.IN.DAT
@@ -3,11 +3,19 @@ i_process_run_mode  = -2
 
 *---------------Constraint Equations---------------*
 * Define number of equality constraints
-n_equality_constraints = 2
+n_equality_constraints = 6
 
 * Equalities
 icc = 1 * Beta
 icc = 2 * Global power balance
+* Tritium particle balance
+icc = 93
+* Deuterium particle balance
+icc = 94
+* Alpha particle balance
+icc = 96
+* Fuelling composition consistency
+icc = 97
 
 * Inequalities
 icc = 5 * Density upper limit
@@ -40,6 +48,26 @@ ixc = 4 * temp_plasma_electron_vol_avg_kev
 boundu(4) = 100.0
 ixc = 6 * nd_plasma_electrons_vol_avg
 
+* 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
+
+
+
+
+
 *---------------Cs Fatigue Variables---------------*
 
 
@@ -310,6 +338,9 @@ temp_plasma_electron_vol_avg_kev       = 12.221383528378944 * volume averaged el
 triang   = 0.5 * plasma separatrix triangularity (calculated if `i_plasma_geometry = 1; 3-5 or 7`)
 i_ind_plasma_internal_norm = 1 * Normalised plasma internal induction selection switch (0 = user input) 
 
+eta_plasma_fuelling = 0.7
+
+
 *----------------------Power-----------------------*
 
 

From 76b2a3ff35d45d1bf5c5af84eace6938cc08f26e Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 3 Aug 2026 10:23:18 +0100
Subject: [PATCH 36/39] Fix bug with nore removed `ovvaref`

---
 process/models/physics/fuelling.py | 6 +++---
 process/models/physics/physics.py  | 1 -
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 3fe44fad05..59fa799a79 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -794,21 +794,21 @@ def output_fuelling_info(self):
             "OP ",
         )
         po.oblnkl(self.outfile)
-        po.ovarrf(
+        po.ovarre(
             self.outfile,
             "Total fuel burn-up fraction",
             "(f_plasma_fuel_burnup)",
             self.data.physics.f_plasma_fuel_burnup,
             "OP ",
         )
-        po.ovarrf(
+        po.ovarre(
             self.outfile,
             "Tritium burn-up fraction",
             "(f_plasma_tritium_burnup)",
             self.data.physics.f_plasma_tritium_burnup,
             "OP ",
         )
-        po.ovarrf(
+        po.ovarre(
             self.outfile,
             "Deuterium burn-up fraction",
             "(f_plasma_deuterium_burnup)",
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 16a4d0054b..4861038169 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -1498,7 +1498,6 @@ def plasma_composition(self):
         # ======================================================================
 
     @staticmethod
-    @nb.njit(cache=True)
     def phyaux(
         fusden_alpha_total: float,
         nd_plasma_alphas_thermal_vol_avg: float,

From dcb960faa41b275b13ea397a1c0863038c20f1d5 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 24 Aug 2026 10:38:54 +0100
Subject: [PATCH 37/39] Fix for new pre-commit errors

---
 process/models/physics/fuelling.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py
index 59fa799a79..78dd95111d 100644
--- a/process/models/physics/fuelling.py
+++ b/process/models/physics/fuelling.py
@@ -17,6 +17,7 @@ def __init__(self):
         self.mfile = constants.MFILE
 
     def run(self):
+        """Run the plasma fuelling model calculations."""
         self.data.physics.molflow_plasma_fuelling_vv_injected_moles = (
             self.data.physics.molflow_plasma_fuelling_vv_injected
             / constants.AVOGADRO_NUMBER
@@ -53,7 +54,7 @@ def run(self):
         )
 
     def output(self):
-        """This model doesn't output to the output file, but it does generate contour
+        """Model doesn't output to the output file, but it does generate contour
         plots of plasma fuel flow rates vs recycling and fuelling efficiency.
         """
 

From 61912451c56eebc13b9b07f080a90983e6a38196 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 24 Aug 2026 13:19:38 +0100
Subject: [PATCH 38/39] Fix expected values in phyaux unit test parameters

---
 tests/unit/models/physics/test_physics.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index cc578d7cfe..15c687ace1 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1850,18 +1850,20 @@ class PhyauxParam(NamedTuple):
     "phyauxparam",
     [
         PhyauxParam(
+            nd_plasma_alphas_thermal_vol_avg=7.499815489237494e18,
             fusden_alpha_total=1.973996644759543e17,
             t_energy_confinement=3.401323521525641,
-            expected_t_alpha_confinement=37.993985551650177,
+            expected_t_alpha_confinement=37.993050845084205,
         ),
         PhyauxParam(
+            nd_plasma_alphas_thermal_vol_avg=7.499880919696304e18,
             fusden_alpha_total=1.9731194318497056e17,
             t_energy_confinement=3.402116961408892,
-            expected_t_alpha_confinement=38.010876984618747,
+            expected_t_alpha_confinement=38.01027347171541,
         ),
     ],
 )
-def test_phyaux(phyauxparam, monkeypatch, physics):
+def test_phyaux(phyauxparam, physics):
     """
     Automatically generated Unit Test for phyaux.
 

From d540f7d9ee584b3e68d0e24d7e36ac374e659170 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 24 Aug 2026 14:05:39 +0100
Subject: [PATCH 39/39] :fire: Remove `burnup_in` variable and related comments
 from physics variables and input files

---
 process/core/input.py                             | 1 -
 process/data_structure/physics_variables.py       | 3 ---
 tests/regression/input_files/st_regression.IN.DAT | 4 ----
 3 files changed, 8 deletions(-)

diff --git a/process/core/input.py b/process/core/input.py
index 0008d256bc..7f8c515546 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -169,7 +169,6 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "beta_vol_avg_min": InputVariable("physics", float, range=(0.0, 1.0)),
     "betbm0": InputVariable("physics", float, range=(0.0, 10.0)),
     "b_plasma_toroidal_on_axis": InputVariable("physics", float, range=(0.0, 30.0)),
-    "burnup_in": InputVariable("physics", float, range=(0.0, 1.0)),
     "radius_plasma_core_norm": InputVariable(
         "impurity_radiation", float, range=(0.0, 1.0)
     ),
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 34251a3946..a23096476f 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -652,9 +652,6 @@ class PhysicsData:
     f_plasma_deuterium_burnup: float = 0.0
     """Deuterium burnup fraction in plasma"""
 
-    burnup_in: float = 0.0
-    """fractional plasma burnup user input"""
-
     b_plasma_vertical_required: float = 0.0
     """Vertical field needed for plasma equilibrium (Bᵥ) [T]"""
 
diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT
index 9f6e8aeb25..f8f90b3f2a 100644
--- a/tests/regression/input_files/st_regression.IN.DAT
+++ b/tests/regression/input_files/st_regression.IN.DAT
@@ -625,10 +625,6 @@ f_p_plasma_separatrix_rad_max  = 0.64
 * DESCRIPTION:   Fraction of alpha power deposited to plasma
 * JUSTIFICATION: Using default = 0.95
 
-*burnup_in = 
-* DESCRIPTION:   User input plasma fuel burnup fraction
-* JUSTIFICATION: Turned off, allow to be calculated in-situ
-
 ixc = 135
 f_nd_impurity_electrons(13) = 0.0005
 *boundl(135) = 0.000495