Particle balance constraint

process/core/input.py

@@ -644,6 +644,15 @@ def __post_init__(self):
     "f_volflow_vac_pumps_impedance": InputVariable("vacuum", float, range=(1e-06, 1.0)),
     "volflow_vac_pumps_max": InputVariable("vacuum", float, range=(1e-06, 1000.0)),
     "molflow_vac_pumps": InputVariable("vacuum", float, range=(0.0, 1e30)),
+    "f_plasma_particles_lcfs_recycled": InputVariable(
+        data_structure.physics_variables, float, range=(0.0, 1.0)
+    ),
+    "eta_plasma_fuelling": InputVariable(
+        data_structure.physics_variables, float, range=(0.0, 1.0)
+    ),
+    "molflow_plasma_fuelling_vv_injected": InputVariable(
+        data_structure.physics_variables, float, range=(0.0, 1e24)
+    ),
     "pflux_plant_floor_electric": InputVariable(
         "heat_transport", float, range=(0.0, 1000.0)
     ),

process/core/solver/constraints.py

@@ -1777,6 +1777,79 @@ def constraint_equation_92(constraint_registration, data):
     )
 
 
+@ConstraintManager.register_constraint(93, "particles/s", "=")
+def constraint_equation_93():
+    """
+    Particle balance
+    """
+    # pscaling: total transport power per volume (MW/m3)
+    # NEED TO ADD PROPER FUSION RATES FOR EACH REACTION
+
+    # Numerator shall be the tritium particle balance
+    numerator = (
+        data_structure.physics_variables.eta_plasma_fuelling
+        * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected
+        - data_structure.physics_variables.fusrat_total
+    ) - (
+        (
+            data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg
+            * data_structure.physics_variables.vol_plasma
+            * data_structure.physics_variables.f_plasma_fuel_tritium
+        )
+        / (
+            data_structure.physics_variables.t_energy_confinement
+            / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled)
+        )
+    )
+
+    # Alpha particle balance
+    denominator = data_structure.physics_variables.fusrat_total - (
+        data_structure.physics_variables.nd_plasma_alphas_vol_avg
+        * data_structure.physics_variables.vol_plasma
+    ) / (data_structure.physics_variables.t_alpha_confinement)
+
+    cc = 1.0 - numerator / numerator
+
+    return ConstraintResult(cc, denominator * (1.0 - cc), denominator * cc)
+
+
+@ConstraintManager.register_constraint(93, "particles/s", "=")
+def constraint_equation_93():
+    """
+    Particle balance
+    """
+    # pscaling: total transport power per volume (MW/m3)
+    # NEED TO ADD PROPER FUSION RATES FOR EACH REACTION
+
+    # Numerator shall be the tritium particle balance
+    numerator = (
+        data_structure.physics_variables.eta_plasma_fuelling
+        * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected
+        - data_structure.physics_variables.fusrat_total
+    ) - (
+        (
+            data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg
+            * data_structure.physics_variables.vol_plasma
+            * data_structure.physics_variables.f_plasma_fuel_tritium
+        )
+        / (
+            data_structure.physics_variables.t_energy_confinement
+            / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled)
+        )
+    )
+
+    # Alpha particle balance
+    denominator = data_structure.physics_variables.fusrat_total - (
+        data_structure.physics_variables.nd_plasma_alphas_vol_avg
+        * data_structure.physics_variables.vol_plasma
+    ) / (data_structure.physics_variables.t_alpha_confinement)
+
+    cc = 1.0 - numerator / numerator
+
+    return ConstraintResult(cc, denominator * (1.0 - cc), denominator * cc)
+
+
+
 def constraint_eqns(m: int, ieqn: int, data: DataStructure):
     """Evaluates the constraints given the current state of PROCESS.
 

process/core/solver/iteration_variables.py

@@ -229,6 +229,18 @@ class IterationVariable:
     174: IterationVariable("triang", "physics", 0.00, 1.00),
     175: IterationVariable("kappa", "physics", 0.00, 10.00),
     176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30),
+    177: IterationVariable(
+        "f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0
+    ),
+    178: IterationVariable(
+        "eta_plasma_fuelling", "physics", 0.01, 1.0
+    ),
+    179: IterationVariable(
+        "molflow_plasma_fuelling_vv_injected",
+        "physics",
+        0.0,
+        1e24,
+    ),
 }
 
 

process/data_structure/numerics.py

@@ -99,10 +99,105 @@ def description(self):
         return self._description_
 
 
-ipnvars: int = 177
+class PROCESSRunMode(IntEnum):
+    """Enumeration of the available PROCESS run modes, which determine the behaviour
+    of the code in various places. This is controlled by the `ioptimz` variable
+    """
+
+    EVALUATION = (-2, "Evaluation mode (no optimisation)")
+    """In this mode, the code will not perform any optimisation, and will instead
+    simply evaluate the constraints for the given input parameters, which is useful
+    for testing and for evaluating the performance of a given design point without
+    trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
+    root-finding method from `scipy.optimize`) to seek a *consistent* solution by
+    varying a subset of the iteration variables until the consistency constraints
+    (equality constraints whose residuals must be driven to zero) are simultaneously
+    satisfied; no figure-of-merit is optimised, and the solver simply tries to find
+    a root of the constraint-residual vector.
+    """
+    OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)")
+    """In this mode, the code will perform optimisation using the VMCON solver
+    (or a custom solver if specified) to try to find a design point that optimises
+    the figure of merit while satisfying the constraints.  This is the default mode
+    of operation for PROCESS.
+    """
+
+    def __new__(cls, value: int, description: str):
+        """Create a new PROCESSRunMode enum member with description.
+
+        Args:
+            value: The integer value of the enum member.
+            description: The description for this run mode.
+
+        Returns
+        -------
+            The new enum member with attached description.
+        """
+        obj = int.__new__(cls, value)
+        obj._value_ = value
+        obj._description_ = description
+        return obj
+
+    @DynamicClassAttribute
+    def description(self):
+        """Return the description for this run mode."""
+        return self._description_
+
+
+class FiguresOfMerit(IntEnum):
+    """Enumeration of the available figures of merit (FoM) that can be used as
+    objective functions for optimisation in PROCESS.
+    """
+
+    MAJOR_RADIUS = (1, "Plasma major radius (R₀)")
+    NEUTRON_WALL_LOAD = (3, "Neutron wall load")
+    P_TF_PLUS_P_PF = (4, "TF & PF coil power")
+    FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)")
+    COST_OF_ELECTRICITY = (6, "Cost of electricity")
+    CAPITAL_COST = (7, "Plant capital cost")
+    ASPECT_RATIO = (8, "Plasma aspect ratio")
+    DIVERTOR_HEAT_LOAD = (9, "Divertor heat load")
+    TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)")
+    TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)")
+    PULSE_LENGTH = (14, "Pulse length")
+    PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor")
+    MIN_R0_MAX_TAU_BURN = (
+        16,
+        "Linear combination of major radius (minimised) and pulse length (maximised)",
+    )
+    NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output")
+    NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit")
+    MAX_Q_MAX_T_PLANT_PULSE_BURN = (
+        19,
+        "Linear combination of big Q and pulse length (maximised)",
+    )
+
+    def __new__(cls, value: int, description: str):
+        """Create a new FiguresOfMerit enum member with description.
+
+        Args:
+            value: The integer value of the enum member.
+            description: The description for this figure of merit.
+
+        Returns
+        -------
+            The new enum member with attached description.
+        """
+        obj = int.__new__(cls, value)
+        obj._value_ = value
+        obj._description_ = description
+        return obj
+
+    @DynamicClassAttribute
+    def description(self):
+        """Return the description for this figure of merit."""
+        return self._description_
+
+
+ipnvars: int = 180
 """total number of variables available for iteration"""
 
-ipeqns: int = 92
+ipeqns: int = 93
 """number of constraint equations available"""
 
 ipnfoms: int = 19

process/data_structure/physics_variables.py

@@ -881,6 +881,14 @@ class PhysicsData:
     molflow_plasma_fuelling_required: float = 0.0
     """plasma fuelling rate (nucleus-pairs/s)"""
 
+f_plasma_particles_lcfs_recycled: float = None
+"""fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)"""
+
+eta_plasma_fuelling: float = None
+"""fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)"""
+
+molflow_plasma_fuelling_vv_injected: float = None
+"""plasma fuelling rate into the vacuum vessel (nucleus-pairs/s)"""
     tauratio: float = 1.0
     """tauratio /1.0/ : ratio of He and pellet particle confinement times"""