Feature/veqpy

documentation/source/fusion-devices/spherical-tokamak.md

@@ -44,7 +44,7 @@ Switch `itart` provides overall control of the ST switches within the code, and
 | --- | --- | --- |
 | `i_plasma_geometry` | 0, 2, 3, 4 | 1, 5, 6, 7, 8 |
 | `i_bootstrap_current` | 1, 2, 3 | 2, 3 |
-| `i_plasma_current` | 1, 3, 4, 5, 6, 7 | 2, 9 |
+| `i_plasma_current` | 0, 1, 3, 4, 5, 6, 7 | 0, 2, 9 |
 | `itfsup` | 0, 1 | 0 |
 
 <center>Table 1: <i> Summary of the switch values in 'PROCESS' that relate to conventional aspect ratio and low aspect ratio machines.</i></center>

documentation/source/physics-models/plasma_confinement.md

@@ -665,6 +665,18 @@ $$
 \tau_{\text{E}} =  0.0821 I_{\text{p}}^{1.02} B_{\text{T}}^{0.11} P_{\text{L}}^{-0.91} \overline{n}_{19}^{0.51} 
 $$
 
+-------------------------
+
+#### 53: Kurskiev spherical tokamak scaling | `kurskiev_st_confinement_time()`
+
+Is selected with `i_confinement_time = 53` [^kurskiev_st]
+
+The fit includes both L-mode and H-mode spherical tokamak discharges and has no isotope-mass term.
+
+$$
+\tau_{\text{E}} =  0.066 I_{\text{p}}^{0.53} B_{\text{T}}^{1.05} P_{\text{L}}^{-0.58} \overline{n}_{19}^{0.65} R^{2.66} \kappa^{0.78}
+$$
+
 -------------------------
 ### Transport Powers
 
@@ -750,6 +762,7 @@ The value of `f_t_alpha_energy_confinement_min` can be set to the desired minimu
 [^22]: G. Verdoolaege et al., “The updated ITPA global H-mode confinement database: description and analysis,” Nuclear Fusion, vol. 61, no. 7, pp. 076006-076006, Jan. 2021, doi: https://doi.org/10.1088/1741-4326/abdb91.
 [^23]: Y. Chen, X. C. Chen, X. F. Wu, and S. Q. Liu, “Energy confinement scaling in the NCST spherical tokamak,” AIP Advances, vol. 16, no. 3, pp. 035043-035043, Mar. 2026, doi: https://doi.org/10.1063/5.0311657.
 [^paz_soldan_neg]: P. Lunia, A.O. Nelson, and C. Paz-Soldan, "Energy Confinement Time Scaling Law Derived from Paz-Soldan NF 2024", doi: https://arxiv.org/abs/2509.04279v2
+[^kurskiev_st]: G. S. Kurskiev et al., “Energy confinement in the spherical tokamak Globus-M2 with toroidal magnetic field reaching 0.8 T,” Nuclear Fusion, vol. 62, no. 1, p. 016011, 2022, doi: https://doi.org/10.1088/1741-4326/ac38c9. The same expression is used for proton-boron spherical tori by H.-S. Xie et al., “ENN's roadmap for proton-boron fusion based on spherical torus,” Physics of Plasmas, vol. 31, no. 6, p. 062507, 2024, doi: https://doi.org/10.1063/5.0199112.
 [^24]: H. Lux, R. Kemp, E. Fable, and R. Wenninger, “Radiation and confinement in 0D fusion systems codes,” Plasma Physics and Controlled Fusion, vol. 58, no. 7, pp. 075001–075001, May 2016, doi: https://doi.org/10.1088/0741-3335/58/7/075001.
 [^25]: H. Lux, R. Kemp, D. J. Ward, and M. Sertoli, “Impurity radiation in DEMO systems modelling,” Fusion Engineering and Design, vol. 101, pp. 42–51, Dec. 2015, doi: https://doi.org/10.1016/j.fusengdes.2015.10.002.
 ‌

documentation/source/physics-models/plasma_current/plasma_current.md

@@ -105,6 +105,15 @@ the switch `i_plasma_current`, as follows:
 
 ---------------
 
+#### User input
+
+Switch value: `i_plasma_current = 0`
+
+The plasma current is taken directly from the input variable `plasma_current_user_input` [A].
+No scaling from $q_{95}$ is applied.
+
+---------------
+
 ### 1. Calculate plasma current shaping function $f_q$
 
 ------------

documentation/source/physics-models/plasma_geometry.md

@@ -269,6 +269,36 @@ $$
 
 ---------------------------------------------------------------------
 
+- `i_plasma_geometry = 13` -- Use a user-provided last closed flux surface (LCFS)
+  as discrete `(R, Z)` points. The following must be input:
+
+  - `r_array` : LCFS radial coordinates (m)
+  - `z_array` : LCFS vertical coordinates (m)
+
+  Arrays may be given as a comma-separated list, e.g.
+
+  ```text
+  i_plasma_geometry = 13
+  r_array = 5.0, 6.0, 7.0, 6.0, 5.0
+  z_array = 0.0, 2.0, 0.0, -2.0, 0.0
+  ```
+
+  From these points PROCESS calculates:
+
+  - `rmajor`, `rminor`, `aspect`
+  - separatrix elongation `kappa` and triangularity `triang`
+    (and the corresponding 95% values via the IPDG89 factors)
+  - poloidal perimeter, surface area, cross-section area and volume by direct
+    contour integration (`cal_integral_geometry`)
+
+  In this mode, `kappa` / `triang` should **not** be treated as independent
+  inputs; they are derived from the LCFS arrays. The derived elongation uses
+  $\kappa = Z_{\max}/a$ with $Z_{\max}=\max(|Z|)$, so for an up-down asymmetric
+  LCFS this is the larger of the upper/lower elongations. Triangularity uses
+  the $R$ coordinate at that same $|Z|_{\max}$ point.
+
+---------------------------------------------------------------------
+
 ### Plasma-Wall Gap
 
 The region directly outside the last closed flux surface of the core plasma is

process/core/constants.py

@@ -125,6 +125,18 @@
 https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Be
 """
 
+M_BORON11_AMU = 11.00930536
+"""Boron-11 atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=B
+"""
+
+BORON11_MASS = M_BORON11_AMU * ATOMIC_MASS_UNIT - 5.0 * ELECTRON_MASS
+"""Boron-11 nuclear mass [kg]
+Atomic mass minus five electron masses.
+Reference: National Institute of Standards and Technology (NIST)
+"""
+
 M_CARBON_AMU = 12.0096
 """Average Carbon atom mass [amu]
 Reference: National Institute of Standards and Technology (NIST)
@@ -232,6 +244,14 @@
 Multiply by the speed of light squared to get the energy released
 """
 
+PB_ENERGY = (
+    (PROTON_MASS + BORON11_MASS) - (ALPHA_MASS * 3.0)
+) * SPEED_LIGHT**2
+"""Proton - Boron-11 reaction energy [J]
+Find the mass difference in the reactants and products of the p-B11 reaction
+Multiply by the speed of light squared to get the energy released
+"""
+
 DT_NEUTRON_ENERGY_FRACTION = ALPHA_MASS / (NEUTRON_MASS + ALPHA_MASS)
 """Deuterium - Tritium reaction energy fraction carried by neutron
 Assuming centre of mass frame as the momenta of the fusion products exceed

process/core/data_structure/base.py

@@ -39,6 +39,7 @@
 from process.data_structure.tfcoil_variables import TFData
 from process.data_structure.times_variables import TimesData
 from process.data_structure.vacuum_variables import VacuumData
+from process.data_structure.veqpy_runtime import VeqpyRuntime
 from process.data_structure.water_usage_variables import WaterUseData
 
 initialise_later = object()
@@ -90,3 +91,4 @@ def __post_init__(self):
         for f in fields(self):
             if getattr(self, f.name) is initialise_later:
                 setattr(self, f.name, f.type())
+        self.veqpy = VeqpyRuntime()

process/core/data_structure/variable_metadata.py

@@ -160,6 +160,11 @@ class VariableMetadata:
     "triang": VariableMetadata(
         latex=r"$\delta_\mathrm{sep}$", description="Triangularity", units=""
     ),
+    "n_lcfs_points": VariableMetadata(
+        latex=r"$N_{\mathrm{LCFS}}$",
+        description="Number of LCFS (R, Z) points",
+        units="",
+    ),
     "f_a_tf_coil_inboard_steel": VariableMetadata(
         latex=r"f_\mathrm{steel}^\mathrm{TF}", description="TF steel fraction", units=""
     ),
@@ -369,6 +374,11 @@ class VariableMetadata:
     "f_nd_impurity_electrons(13)": VariableMetadata(
         latex=r"$Xe_{\mathrm{f}}$", description="Impurity fraction (Xenon)", units=""
     ),
+    "f_nd_impurity_electrons(15)": VariableMetadata(
+        latex=r"$B_{\mathrm{f}}$",
+        description="Boron fuel fraction (when proton-boron fusion)",
+        units="",
+    ),
     "pdivmax_over_rmajor": VariableMetadata(
         latex=r"$P_{\mathrm{div}}/R_\mathrm{maj}$ [MW/m]",
         description="Divertor power per major radius",

process/core/init.py

@@ -34,6 +34,7 @@
 from process.data_structure.stellarator_variables import StellaratorModel
 from process.data_structure.superconducting_tf_coil_variables import TFWPIntegerTurnType
 from process.models.pfcoil import PFLocationTypes
+from process.models.physics.plasma_current import PlasmaCurrentModel
 from process.models.physics.profiles import (
     DensityProfilePedestalType,
     PlasmaProfileShapeType,
@@ -376,20 +377,53 @@ def check_process(inputs, data):  # noqa: ARG001
         )
 
     #  Fuel ion fractions must add up to 1.0
-    if (
+    if data.physics.i_fusion_reactions == "p-b11":
+        if (
+            abs(
+                1.0
+                - data.physics.f_plasma_fuel_boron11
+                - data.physics.f_plasma_fuel_proton
+            )
+            > 1e-6
+        ):
+            raise ProcessValidationError(
+                "p-b11 fuel ion fractions do not sum to 1.0",
+                f_plasma_fuel_boron11=data.physics.f_plasma_fuel_boron11,
+                f_plasma_fuel_proton=data.physics.f_plasma_fuel_proton,
+            )
+        if (
+            abs(
+                data.physics.f_plasma_fuel_deuterium
+                + data.physics.f_plasma_fuel_tritium
+                + data.physics.f_plasma_fuel_helium3
+            )
+             > 1e-6
+        ):
+            raise ProcessValidationError(
+                "D-T Fuel ion fractions do not sum to 0.0",
+                f_plasma_fuel_deuterium=data.physics.f_plasma_fuel_deuterium,
+                f_plasma_fuel_tritium=data.physics.f_plasma_fuel_tritium,
+                f_plasma_fuel_helium3=data.physics.f_plasma_fuel_helium3,
+            )
+    elif (
         abs(
             1.0
             - data.physics.f_plasma_fuel_deuterium
             - data.physics.f_plasma_fuel_tritium
             - data.physics.f_plasma_fuel_helium3
         )
-        > 1e-6
+        > 1e-6 or abs(
+            data.physics.f_plasma_fuel_boron11
+            + data.physics.f_plasma_fuel_proton
+        ) > 1e-6
     ):
         raise ProcessValidationError(
             "Fuel ion fractions do not sum to 1.0",
             f_plasma_fuel_deuterium=data.physics.f_plasma_fuel_deuterium,
             f_plasma_fuel_tritium=data.physics.f_plasma_fuel_tritium,
             f_plasma_fuel_helium3=data.physics.f_plasma_fuel_helium3,
+            f_plasma_fuel_boron11=data.physics.f_plasma_fuel_boron11,
+            f_plasma_fuel_proton=data.physics.f_plasma_fuel_proton,
         )
 
     if data.physics.f_plasma_fuel_tritium < 1.0e-3:  # tritium fraction is negligible
@@ -403,6 +437,15 @@ def check_process(inputs, data):  # noqa: ARG001
             stacklevel=2,
         )
 
+    if (
+        data.physics.i_plasma_current == PlasmaCurrentModel.USER_INPUT
+        and data.physics.plasma_current_user_input <= 0.0
+    ):
+        raise ProcessValidationError(
+            "plasma_current_user_input must be positive when i_plasma_current=0",
+            plasma_current_user_input=data.physics.plasma_current_user_input,
+        )
+
     if data.impurity_radiation.f_nd_impurity_electrons[1] != 0.1:  # noqa: RUF069
         raise ProcessValidationError(
             "The thermal alpha/electron density ratio should be controlled using"
@@ -495,6 +538,11 @@ def check_process(inputs, data):  # noqa: ARG001
                 data.numerics.boundu[3], data.numerics.boundl[3]
             )
 
+        if data.physics.i_equilibrium_solve == 1 and data.physics.i_alphaj != 0:
+            raise ProcessValidationError(
+                "i_equilibrium_solve=1 requires i_alphaj=0 (user alphaj for veqpy j_tor)"
+            )
+
         # Density checks
         # Issue #589: Pedestal density is lower than separatrix density
         pedestal_type = DensityProfilePedestalType(
@@ -529,6 +577,21 @@ def check_process(inputs, data):  # noqa: ARG001
                     }
                 ),
             )
+        if (
+            pedestal_type == DensityProfilePedestalType.GREENWALD_FRACTION 
+            and data.physics.i_equilibrium_solve == 1
+        ):
+            raise ProcessValidationError(
+                "Pedestal and separatrix densities must be input as absolute "
+                "values (nd_plasma_pedestal_electron, "
+                "nd_plasma_separatrix_electron) with "
+                "i_nd_plasma_pedestal_separatrix = 0 when i_equilibrium_solve = 1. "
+                "Greenwald-fraction inputs are not allowed.",
+                i_nd_plasma_pedestal_separatrix=(
+                    data.physics.i_nd_plasma_pedestal_separatrix
+                ),
+            )
+
 
         if (
             abs(data.physics.radius_plasma_pedestal_density_norm - 1.0) <= 1e-7
@@ -641,6 +704,14 @@ def check_process(inputs, data):  # noqa: ARG001
                 "and is recommended for the Reinke model",
                 stacklevel=2,
             )
+    if (data.physics.i_fusion_reactions == "p-b11") and (
+        data.physics.i_l_h_threshold in {1, 2, 3, 4, 5, 15, 16, 17, 18}
+    ):
+        logger.warning(
+            "p-B11: i_l_h_threshold has no z_eff or m_ion correction; "
+            "a scaling that includes either is recommended",
+            stacklevel=2,
+        )
     i_single_null = DivertorNumberModels(data.physics.i_single_null)
     match i_single_null:
         case DivertorNumberModels.DOUBLE_NULL:

process/core/input.py

@@ -17,7 +17,10 @@
 from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.data_structure.numerics import N_ITERATION_VARIABLES_MAX
 from process.data_structure.pfcoil_variables import N_PF_GROUPS_MAX
-from process.data_structure.physics_variables import N_CONFINEMENT_SCALINGS
+from process.data_structure.physics_variables import (
+    N_CONFINEMENT_SCALINGS,
+    N_LCFS_POINTS_MAX,
+)
 from process.data_structure.scan_variables import IPNSCNS, IPNSCNV
 
 if TYPE_CHECKING:
@@ -48,6 +51,28 @@ def _icc_additional_actions(
     data.numerics.n_constraints += 1
 
 
+def _lcfs_array_additional_actions(
+    _name, value, array_index, _config, data: DataStructure
+):
+    """Track how many LCFS (R, Z) points were provided in the input file.
+
+    Raises
+    ------
+    ProcessValidationError
+        If an LCFS array exceeds the maximum allowed length.
+    """
+    if isinstance(value, list):
+        n_points = len(value)
+        if n_points > N_LCFS_POINTS_MAX:
+            raise ProcessValidationError(
+                f"LCFS array length {n_points} exceeds maximum "
+                f"{N_LCFS_POINTS_MAX} ({_name})"
+            )
+        data.physics.n_lcfs_points = n_points
+    elif array_index is not None:
+        data.physics.n_lcfs_points = max(data.physics.n_lcfs_points, int(array_index))
+
+
 @dataclass(slots=True)
 class InputVariable:
     """A variable to be parsed from the input file."""
@@ -435,6 +460,9 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "f_rad": InputVariable("stellarator", float, range=(0.0, 1.0)),
     "f_sync_reflect": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_plasma_fuel_tritium": InputVariable("physics", float, range=(0.0, 1.0)),
+    "f_plasma_fuel_boron11": InputVariable("physics", float, range=(0.0, 1.0)),
+    "f_plasma_fuel_proton": InputVariable("physics", float, range=(0.0, 1.0)),
+    "f_nd_protons_electrons_input": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_beam_tritium": InputVariable("current_drive", float, range=(0.0, 1.0)),
     "f_vforce_inboard": InputVariable("tfcoil", float, range=(0.0, 1.0)),
     "f_w": InputVariable("stellarator", float, range=(0.1, 1.0)),
@@ -1005,16 +1033,24 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "i_density_limit": InputVariable("physics", int, range=(1, 8)),
     "i_diamagnetic_current": InputVariable("physics", int, choices=[0, 1, 2]),
     "i_div_heat_load": InputVariable("divertor", int, choices=[0, 1, 2]),
-    "i_l_h_threshold": InputVariable("physics", int, range=(1, 21)),
+    "i_l_h_threshold": InputVariable("physics", int, range=(1, 24)),
     "i_pf_current": InputVariable("pf_coil", int, choices=[0, 1, 2]),
     "i_pfirsch_schluter_current": InputVariable("physics", int, choices=[0, 1]),
-    "i_plasma_current": InputVariable("physics", int, range=(1, 9)),
-    "i_plasma_geometry": InputVariable("physics", int, range=(0, 12)),
+    "i_plasma_current": InputVariable("physics", int, range=(0, 9)),
+    "plasma_current_user_input": InputVariable("physics", float, range=(1.0e5, 1.0e8)),
+    "i_plasma_geometry": InputVariable("physics", int, range=(0, 13)),
     "i_plasma_shape": InputVariable("physics", int, choices=[0, 1]),
     "i_plasma_wall_gap": InputVariable("physics", int, choices=[0, 1]),
+    "r_array": InputVariable(
+        "physics", float, array=True, additional_actions=_lcfs_array_additional_actions
+    ),
+    "z_array": InputVariable(
+        "physics", float, array=True, additional_actions=_lcfs_array_additional_actions
+    ),
     "i_pulsed_plant": InputVariable("pulse", int, choices=[0, 1]),
     "i_q95_fixed": InputVariable("constraints", int, choices=[0, 1]),
     "i_r_cp_top": InputVariable("build", int, choices=[0, 1, 2]),
+    "i_calculate_radiation": InputVariable("physics", int, choices=[1, 2]),
     "i_rad_loss": InputVariable("physics", int, choices=[0, 1, 2]),
     "i_shield_mat": InputVariable("fwbs", int, choices=[0, 1]),
     "i_single_null": InputVariable("physics", int, choices=[0, 1]),
@@ -1042,6 +1078,8 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "ifetyp": InputVariable("ife", int, range=(0, 4)),
     "ifueltyp": InputVariable("costs", int, choices=[0, 1, 2]),
     "i_plasma_ignited": InputVariable("physics", int, choices=[0, 1]),
+    "i_fusion_reactions": InputVariable("physics", str, choices=["d-t-he3", "p-b11"]),
+    "i_nd_plasma_protons": InputVariable("physics", int, choices=[0, 1]),
     "i_blkt_module_segmentation": InputVariable("fwbs", int, choices=[0, 1]),
     "inuclear": InputVariable("fwbs", int, choices=[0, 1]),
     "iohcl": InputVariable("build", int, choices=[0, 1]),
@@ -1053,6 +1091,7 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "i_beta_norm_max": InputVariable("physics", int, range=(0, 5)),
     "i_ind_plasma_internal_norm": InputVariable("physics", int, range=(0, 2)),
     "i_alphaj": InputVariable("physics", int, range=(0, 1)),
+    "i_equilibrium_solve": InputVariable("physics", int, choices=[0, 1]),
     "i_fw_blkt_shared_coolant": InputVariable("fwbs", int, choices=[0, 1, 2]),
     "ireactor": InputVariable("costs", int, choices=[0, 1]),
     "irefprop": InputVariable("fwbs", int, choices=[0, 1]),